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1. INTRODUCTION 
A. Overview 

Time-accurate Implicit finite-difference schemes for the Euler nnd com- 
pret!'»lble Navler-Stokes equations are used to obtain steady as well as 
unsleady flow-field solutions. If only a steady-state solution Is required, 
Iterative paths that are not restricted to be time accurate can be sought 
to acceleiutc steady-atate convergence. This Is the concept of relaxation 
which has been used succeesfully for Invlscld transon’c flow. 

The current time accurate Implicit algorithms [1-6] for the Euler or 
compressllite Navler-Stokes equations rely on approximate factorization or 
alternating direction (ADI) techniques to achieve computational efficiency. 
The same technique Is the basis of many of the most successful relaxation 
pro.edurea (e.g., [7-Jl]). As a consequence, It would seem that Implicit 
algorithms developed for time accurate flow simulation could be adapted 
Into succossful relaxation procedures, and indeed this Is the case. 

In tills work, tlu* Iterative convergence properties of a currently 
popular approximate-factorization Implicit finite-difference algorithm are 
studied botli analytically and experimentally. These studies are limited to 
the two-dimensional Euler equations, with emphasis on transonic flow compu- 
tations. However, tho major results are expected to apply to those flows 
that are t'overned by the complete Navler-Stokes equations, but In which the 
convection phenomena atlll plav the most important role In the determination 
of the esuentlal features of the numerical algorithm, at lea^* from the 
standpoint of stabllliv. 
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To achieve better numerical efficiency, large t Ime-ateps are often 
needed for problems that are unduly stiff. To permit this, modifications 
to the algorithm are nade (In Section IB) In an attempt to enhance Its sta- 
bility properties. The success of this attempt is supported by a theoreti- 
cal analysis (Chapter II), and a numerical experimentation (Chapter III). 
With achievement of stable large time-steps permitted, another technique is 
also employed to Improve the iterative convergence rate. This technique, 
which consists of using a cyclic sequence of time-steps, appears promising 
after examination of a simple model problem (see Section IIC). In 
Chapter III, a variety of numerical experiments are conducted on the modi- 
fied algorithm and the use of a sequence of time-steps. 

Finally, it was observed in the course of this work, that the numerical 
algorithm could be subject to a particular form of Instability due to 
variable coefficients. A discussion on this topic Is presented In 
Chai»ter IV. 

In Section IB, which follows, the definition of the base algorithm is 
recalled, and a modified algorithm Is proposed. 


B. Governing Equations and Numerical Algorithm 
The conservative form of the Euler equations In Cartesian coordinates 
and for two-dimensional flow is given by: 

3^q + + 3y? - 0 (1) 


wliere 


q - 



i - 



, and 
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Implicit finite-difference schemeH developed for the Euler equetions in 
nonconservatlve (11 or conservative fo <n (2-3) share the essential features 
of central spacial differencing (for stability) and alternatlng-dlrectlon- 
llkf structure (for efficiency). The conservative differencing scheme Is 
used here for transonic flow applications because It correctly captures 
shock waves, but the results obtained here should apply to the nonconserva- 
tive differencing scheme as well. 

The implicit f Inl t e-dl ffc-encing scheme can be represented as 

(I + h6^")(I + h6yB")(q"‘^‘ - q") 

- -At(5j" + (Syl") - + (VyAyl^lq" (2) 

where 


and Ay are central three-point difference operators 


A, B are the Jacobian matrices 


ft) ■ [*] 


h • I' At, with 0 ■ I or 1/2 for Euler implicit, or trapezoidal time 
dlf ferencinj; 

(VA)'^” ere fourth-order numerical dissipation terms with coefficient 


•e i 


Mr , *n 


q • with Xj • (J - l)Ax and * (k - D^y. 

The operators A and VA are understood to operate on any product of terms 
that follow to their right and, for example. 


q 




+ Mjk - *<ij+,,k + <ij+. _k 


Central difference operators are used because A and B usually have 
botli positive and negative eigenvalues, for one sign at v»hlch the algorltltm 
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Is always unstable for only forward or only backward spaclal differencing. 

The fourth-order numerical dissipation terms provide damping and art needed 
to control what is usually referred to as nonlinear instability. They also 
serve to smooth out small numerical inconsistencies especially due to 
slightly improper boundary conditions, 

Although the basic differencing is stable for linear problems without 
dissipation added, the scheme given by Equation (2) has to be modified if 
the dissipation is to be allowed to increase with the use of large values 
of At. More precisely, the dissipation coefficient should vary 

directly as At to prevent nonlinear instability and to maintain steady- 
stale consistency. It is only in this way that the steady-state solution 
can be independent of At. However, because the numerlc.al dissipation is 
added explicitly, use of > 1/16 would In Itself cause linear instabil- 
ity, Consequently, Cj, cannot be maintained proportional to At and very 
large values of At cannot be taken without effectively reducing the amount 
of added numerical dissipation. In many flow-field problems, lack of suffi- 
cient numerical dissipation can cause instability. 

Adding numerical dlsalpatlon Impllcltlv would allow to assume anv 

positive value, and in particular, tg could vary with At. Unfortunately, 
use of fourth-order lupllclt numerical dissipation requires the inversion 
of block pentadlagonal matrices which are twice as costly as the block tri- 
diagonal inversions required for Equation (2), Use of second-order smoothing 
all('W8 tri diagonal structure but is Inaccurate. However one expects that 
lean restricted values of Cg could be obtained if a proper portion of the 
ouiik rical dissipation that fits within the tridlagonal structure is treated 


*»a ‘ 
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Implicitly. Thin concept ultimately leads to the following modification of 
the numerical algorithm 

(I + h6jjA" - + h«yj" - 

- -At (6^1" + 6y?") - 

which has now been Imfilemented in recent flow-field codes [5,6]. 

In the follovdng chapters, this modified algorithm is analyzed and com- 
• pared to the- base dlfierencing scheme. Equation (2), from the atandpolnt of 

Iteiatlve convergence. 


,!ISS! SBWi ns*»^4 
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II. STABILITY ANALYSIS AND APPLICATIONS 

The efficiency of a finite-difference algorithm for tlme-marchlng 
problems depends crucially on the stability limitations this algorithm is 
subject to. In Section TIA, the basic notions related to stability are 
recalled. In Section IIB, stability bounds are derived for the case where 
the numerical algorithm is applied to a scalar, linear, partial-differential 
equfition, with constant coefficients and linear boundary conditions. Appli- 
cations of the results of this analysis are considered in Section IIC. 

A. Generalities 

The importance ol the stability condition is reviewed hero for finite- 
difference- methods in general, and for relaxation techniques in particular. 

Recall that a numerical algorithm, for an initial-value problem, is 
said to be stable when, for arbitrary bounded starting solution u° ■ u(0), 
the solution u” produced by n applications of this algorithm remains 
bounded an n tends to infinity. This limit may result from ronslderlng 
either one of two limiting processes which are: (1) a mesh refinement 

process, and (2) a search for a steady-state solution. 

In the case of a mesh refinement process, one evaluates a sequence of 
solutions u ^ (1 “ 1,2, . . .) which are all candidate approximations to 
the exact solution u(tj) of the initial-value problem, for some fixed 
flniil time tf. At the 1th step in this process, nj^ applications of the 
algorithm are made, with initial solution u° ■ u(0) and using a time-step 
Ati - tf/nj^ which tends to zero as 1 tends to Infinity. P. D. Lax (see 
Rlchtmyer and Morton (121) has shown that given a properly posed Initial- 
value problem and a finite-difference approximation to It that satisfies 
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the crtnttl ► r i»ncy condition, the stability condition Is the necessary and 
Hufllcient condition lor convergence. Here the convergence Is the one oi 

- u(tf)l to zero as 1 tends to Infinity (or At^ tends to zero); 
that is. the convergence ol the finite-difference integral operator to the 
exai t integral operator over a fixed domain in the limit of a mesh 
ref inement , 

In the case of the search for a steady-state solution* the time-step 
\t nuw bt fixed and n tends lo Infinity because the final time • n At 
should do so. There the stabllitv condition Is not sufficient for (steady- 
state) convergence , and one usually relies on numerical evidence to demon- 
strate tht latter. 

For both cases, violation of the stability condition produces an ampli- 
fication i»f the various forms of errors that are present In the numerical 
solution. These are: truncation errors (due to Inexact differentials), 

rouud-off errors (due to truncated arithmetics), errors due to slightly 
Inconsisttnt boundarv conditions, etc. Fox linear (constant coefficient) 
algi'r i thmt. , the growth of the errors, !f it happens, is generally exponen- 
tlvii with n (sometltms polynomial . or a combination of the two), so that 
the luimevical solution verv rapldlv becomes fotallv meaningless whenever 
computable. However, most schemes are stable when operating with a time- 
step that does not exceed a certain maximum allowable value 
unfittunat ely decreases with the mesh spacing parameters Ax and Av, and 
alsi* depends (for nonlinear sclienn's) on the solution u*^ Itself. For 
example, lor usual ex|>llcit algorithms (e.g.. (13)1 applied to the Kuler 
equiitions, stabllitv is enfc cd bv the well-known Courant , Friedrichs, 

^ew^ (CFL) condition ll^]. T. s condition conslderablv reduces the 



8 


efficiency of these algorithms when used as relaxation techniques. This is 
particularly true when, for an accurate resolution, a very fine mesh Is 
required. On their part, implicit algorithms usually have the favorable 
property of being unconditionally stable, at least for some simple test 
equations. In practice, such unconditionality is rarely truly achieved, 
but time-steps that are significantly larger than those permitted by the 
CFL condition can be successfully used (see Chapter III). This Is the 
reason that motivated the choice of an implicit algorithm in this work on 
relaxation . 

These considerations suffice to explain the importance of stability for 
finite-difference time-marching techniques. However, when such a technique 
is (‘raploycd as an artifice to solve a problem where time does not appear, 
one might want to relate the stability condition to the assumption of a 
(known) tlieorem dealing with relaxation techniques per se. This is the 
'‘contraction-mapping theorem” (see, e.g., [15]) which can be stated as 
follows: Given a closed domain D in a complete normed vector space (e.g., 
R®), and an applicaticm f, with domain D and range included in D, which 
Is contracting in the sense that 

Vu,v £ D , if(u) - f(v)il < pHu - vK (4) 

for some real positive number p 1, the following statements are true: 

(1) the equation u ■ f (u) has a unique solution u* (on D), and (2) for 
any € D, the sequence given by and » f(u^) 

(n 0,1,2, . . . ) is well defined and convetges to u’*. Also, the fol- 
lowing bound holds: 

P t n n— 1 ■ 

1 I U - U ■ 

1 * p 


lu" - U*l 


(5) 
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Q 

When a f inir forer.ce inethod ir? used as a lelaxation technique, the 
Iterative formula can indeed be written as 

u ^ f(u ) (6) 

where f(u) can i^en*»rall' be cast Into the following quaei-llnear form: 

f(u) • L(u)u b(u) (7) 

where •* u(n^t) Is the (m-dimenaional ) solutlon-vectcr , L(u) Is an 
m m coefficient matrix and b(u) is an m-vector generally resulting from 
the appliiatton of boundary conditions. T\w enfortunate dependence on u 
of L(u) and b(u) renders the analysis very difficult In very general cases. 
For this reason, one is generally satisfied when successful in proving 
sta!>llity o.' the algorithm in the special case where the co^*f f ic ients are 
frozen to some, perhaps arbitrary but fixed, nominal values L and b. Then, 
the boui^dedness of u^, for arbitrary u^. Is equivalent to the following 
8ta)>lllty condition: 

•U ^ 1 (8) 

♦or some norm. Clear! v, this condition Is a weak form of the assumption 
that f IS contracting of the cited theorem (see Kqiut Ion (4)^. 

If the matrix L can be diagonalized. It Is convenient to use tlie 
spei'tral norm for which Kquatlon (8) reduces to: 

*^(L) i I (9) 

Matrices that are involved In finite-difference equations are usuallv band 
matrices. Despite this simplification, the determination of the eigensystem 
of L Is usually dlfticult. For this reason, most analyses apply to simple 
test cases, such as tite one of a scalar, linear part lal-dlffere.it ial equa- 
tion. with coefficients assumed constant in both time and space, and for 
so»i‘ simple boundary ionditions (generally periodic sometimes fixed, rarely 



t 
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mor:* sophistic&trd) . In Section IIB, the stability analysis is developed 
for this simple case. Such analysis can be invalidated by either on*i of the 
following realities: 

1 . Dimensionality 

2. Nonlinearity or time-dependence 

3. Spacial dependence of the coefficients matrices 

4. Complex boundary condition procedures 

In Chapt'-t IV, a form of instability due to spacial variation of the 
Jacobian matrices of the Euler equations (item 3) is discussed - 

B. The Case of a Scalar Linear Equation 


1 . Generalities 

If the Jacobian matrices A and B of the Euler equations commuted and 
wero constant, the governing equations could be diagonalized into four 
scalar equations of the form: 

u^ + au^ + bUy = 0 (10) 

which is the first-order wave equation in two dimensions. Although these 
hypotheses are not satisfied by the Euler equations, the case of applica- 
tion of the numerical algorithm to Equation (10) is expected to reveal the 
essential properties of this algorithm. For this case, if linear boundary 
conditions^ are assuim^d, it is convenient to rewrite the finite-difference 
equ;ition (Equation (3)) in the following matrix form (derived in Appendix A): 
® - u") - -(B^ ® ly + © By)u" (11) 


^The analysis will be made for periodic or specified boundary 
conditions. 



I J t t I ♦ 


» ^ ^v4 ^ *4 


whtMC the foliowluK ut flnltlonr. b.nve been uned 


'x ^ » ‘ i\ 


\ - ^‘x^'x ^ 

- (2Ax)Sj^ - Ti ld(-l,0,l) 

(12) 

l\ - \\^^ - Trld(-U2,-1) 

Dj^ " ^^x • aecond or fourth-order smi’»othlng 

- aAt/(2A8) , half of a Couranl number ^ 

and Ay% \% ^«y* 0^,, P*,, ami Vy are defined In a similar way. For a mesh 
contalniuK. d'^K Intel lor ^rid points, x-suhsor Ipt ed and y-aubscr Ipted 
mati lees are of dimension ami K^K, respectively. It is assumt'd that 

the J'^K components i»f the solution vector u*' are conventionally ordered 
as lollows: 

n . n n n n n n n n n . t , . ^ v 

u • ‘ • • 

where as usual • u(Xj, yj^, ’^n^* also assumed that this vector 

has been defined In such a way tliat u'' ■ 0 Is the solution of the differ- 
ence' equation t hf t one hopes to attain at the steady stale. The homojteneltv 
of liquation (11) results from this Implicit convention. Flnallv, defini- 
tions Hnd eigensvstems of trldlagonal matrices for varlovis boundary condi- 
tions are given In Appendix B. 

Inveiting Equation (11) yields the new equation: 

n*f I n , , . 

u - Lu tlA) 

.A\ou' I, « 1 - - A~* «f>A~'By. 

Lot X and Y he two nonatniitulHr irntrlooH of alroa .1''.! and K^K 


roHpoct ivo I y , to ho choaen lator. Upon defining 
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v" = (X ® Y)”'u" (15) 

and makinp, various applications of Equations (A3) through (A5) (see 
Appendix A), Equation (14) becomes: 

= Av" (16) 

where the matrix A, which is defined by 


A - (X ® Y)“'l(X ® Y) 

- I - ®y"1a^,'y - X~U~^X ®Y"'A“'R^Y 

- I ® A^‘ - A"1 (17) 


in which, for example « 

Ax - 

&x * x"^B^y 


(18) 


Is slmllai' to the matrix L. 

Observe that for the simple case where no smoothing is applied 
(Ec * ^i matrix A can be reduced to a diagonal form. For this, 

it .suffici s to choose X and Y to diagonal :fze the matrices and C„. 

A > 

Now. considering the general case where and are nonzero, and assum- 
ing the nuitrix A ditigonalizable, an inspection of Equation (14) or (16) 
Indicates that the solution u'^ (alternately v”) is bounded for arbitrary 
starting rolutlon, if and only if the spectral radius of the matrix L 
(alternately A) is less than or equal to unity. I,, is desired to determine 
the conditions on fg, and 0 under which this requirement is met for 
arbitrary values of the parameters and \>^ that control the time-step 
At ( unconditional stiibillty") . This Is done in the next two sections for 


some assumed bounl ry conditions. 




•^ ' I 1 
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2. The cawe of per iodi c boundary c onditions ' 

The ( ase of periodic boundary conditions Is of Interest because It 
permits tlie development of a rlKorous analysis of a pure initial-value 
problem. However, as one can anticipate by observlnj? that the exact solu- 
tion for I his problem Is given by 

u(x,y,t) • u(x - at, y - bt, 0) (19) 

It tioes m't provide a satisfactory test case for studying steady-state con- 
vergence. Neverthelens, the analysis of this case will be performed here 
us .1 guideline for the treatment of another case. 

When periodicity conditions are applied. It Is convenient to assume 

that u” Is defined for all (positive or negative) integer values of 
■Ik 

J and k and that 

Vvd.k+pK " “jk " integers) (20) 

The forward and backward shift operators (acting on either J or k) are then 
Invi'rse oi one another (see Appendix A); ao are their nuUrlx ropreaentat Iona 
whl» h thus, can he slimil tane.^usl v diagonalized. As a consequence, the 
matrices in Equation (11) with the same subscript (x or y) which are llneiir 
cornl* Inat tons of their powers are also diagonalized bv the same transforma- 
tion. The (clrcuiant) eigenvectors of these matrices are then choaen to 
construct X nd Y (for detalla sec Appendix B) . Thia gives; 


Beam [16] originally showed the unconditional stability of the algo- 
ritlim whan applied (consistently) to the e<|uatlon: 

% “x + “v * • <“xx + “yy^ 

rhls analysis Is extended here to the case where the dissipation terms may 
be lourth derivatives as well as second derivatives. Also. In the present 
analvsl.s, these dissipation terms are differenced in a wav that Is not 
nec»*ssarl 1 v t iipe-’accur ^te. 
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■ l/»^ fxp(mlOj) 
- 1//K exp(mlO^' 


Cl) 


where Oj - 2 ti() - 1)/.I (J “ 1.2, . . . ,J) and 0]^ ■ 2it (k - 1)/K 
(k - I,:, . . . ,K).^ The eigenvalue of the J forward shift operator asso- 
ciated to the eigenvector Xj ■ (X^j) la simply ” exp(l('j). 

Similarly, the eigenvalue of the k forward shift operator associated to 
the eigenvector Yj^ - exp(10j^). The eigenvalues of other oper- 

ators are obtained as linear combinations of the powers of exp (10 j) or 
exp(lOj^). For example, the eigenvalues of the matrices A^, C^, 

and are given, respectively, by: 

- 1 + 9Vjj(iCj) + c^dj 

v„(lc.) + e.d. 


M 'x'-j' • “i”j 
ICj ■ exp (10^) - exp(-10j) ■ 21 sin 0^ 

dj " -exp(-jOj) + 2 - exp(iOj) ■ 2(1 - cos 0^) 

dj " ^ ~ second or fourth-order smoothing 


( 22 ) 


The eigenvalues of the matrices Ay, By, Cy, Dy, and Dy are denoted by 
*k’ '’k* ^* k’ ^k’ '*k’ given by an equation similar to Equa- 

tion (22). 


^ Remark : It is shown, in Appendix B, that X and Y are unitary, that Is 

X"^ - y/ (adjoint of X) - X*^ 

Y"' - Y* (adjoint of Y) - Y^ 

This is because X and Y represent (flnlte-dlmenslonal. Inverse) Fourier 
transforms. 




I ' 


I 
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Wiih Liiirt iiioice of X and Y, the matrices A^* and By in 

Kqu;illon (17) are all diagonal. As a result, the matrix A Itself oocoraes 
a diagonal matrix with eigenvalues given by; 

'jk " ‘ ■ ‘“j''’j’*k‘ ■ ‘’r<‘k''>k' 

*i‘k - ‘’>.1 * "k* 


This gives: 

''jk ■•■ 

Sk " a~rr/~ ^24) 

in which and arc real and given by: 

“jk “ ' '■i‘*k^ ' 

t<jk ■ 

I (25) 

“jk ’ “jk “ 'e<‘^] + 

V • "jk - ■*■ 

whei'f C| • VjjCj and c.. ■ stability condition (U^i^l i D then 


becomea : 


(U^ 4 ’ fu^ 4 

Jk Jk - “jk ‘jk 


■^‘jk‘o<‘*j + "o' 4 c.) 4 (c, 4 c.y < 0 

or 

4 20 (i j 4 Cp)[(l 4 ' jdj^)Cj 4 (1 4 i^d^Jc^l - (c, 4 > 3 



1 
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or 

T + <KcpC2) - 0 (2,) 

where 

T • Cg(d^ + d,^)[2(l + idj)(l + e^dj^) - c^CtlJ + d^l (28) 

and Q(cj,C 2 > Is a quadratic form in CpCj given by; 

Q(cj,C 2> - [26(1 + r^dj^) - Ucj^ + 2[0(2 + t^Uj + f^d^) 

- 0^c„(dj + dp - l]c,C2 + [20(1 + e^dp - 1]C2^ (29) 

Note* that If no smoothing Is appll^'d (c^ “ “ 0), T vanishes identically, 

while Q(C|,C 2 ) reduces to (20 - l)(cj + From this, one concludes 

that the I uler explicit method (0 » 0) is unconditionally unstable for this 
cast-, while the trapezoidal time-differencing method (0 * 1/2) as well as 
the Euler implicit method (0 * 1) are both unconditionally stable, as well 
known [2,3]. 

Now consider again the general case where and ^ are nonzero. Since 

the wave speeds a and b as wt*ll as the time-step At ought to be arbitrary, 
the parameters and Vy and consequently Cj and must be considered as 
free parameters. The condition expressed in Equation (27) then breaks into 
two: 


T > 0 

(30) 

Q(c 1 .Cj) > 0 

(31) 


Equation (31) will be examined first. In view of Equation (29) it is appar- 
ent that its satisfaction requires, in particular, that the coefficients of 
Cj^ and Cp^ in Q(cj,C 2 ) be nt»nne^’:ative, Tills gives the following necessary 
conditions : 





% 




I 



i 


1 


I 
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20(1 + > 1 
20(1 + tjdj) i 1 


O.') 


The definitions of dj and dj^ (Equation (22)) Indicates that these eigen- 
values are positive (dissipation), but tend to zero (for fixed values of 
J and k) In the limit of a mesh refinement Ax, A * -► 0 (or J,K ► «')• Hence, 
Equation (32) requires that 



(33) 


Sufficiency is <)btalned by enforcing, also, that Q(cj,c.) be nonfactor- 
able. This gives the following condition: 

(0(? + f^dj + " 0?r^(dj + dj^) - ll^ 

- f20(l + fjd^) - 1](20(1 + r^dj) -1)^0 (34) 

Expanding above quartlc form in 0 would reveal that 0^ can be factored 
In It. For this reason, after a few simpllflcaclons, a condition equivalent 
to I-quatlon (34) can be obtained In the form: 

gj^(fl) i 0 (35) 

where 8 jj 5 (*^) ® quadratic form in 0 given by: 

gjj^(O) - 0-'^e^^(dj + dp2 - 20fg(dj + d^)(r^dj + r^d^ + 2) 

+ 2 fg(dj + dp + L^^(dj - (36) 

Equation (35) roust be enforced for all values of J and k and for the par- 
ticular value of 0 which corresponds to the chosen method (0 ■ 1/2 for 
trapezoidal timu-dlfferenclng, 0-1 for the Euler implicit method). In 
this way, a condition on Cg and results. Before expllcltlng this con- 
dition, assume momentarily that r^, r^, and 0 have been chosen to satisfy 
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precisely this condition. Then, corfalnlv, for some arbitrarily choaen 
values of J and k, It is true that 

if one defines g*j^ - Min achieves its minimum for 

0 • 2)/[t.g(dj + dpi sc that; 

*Jk ■ ^ f 2e,(<ij + d’) H- - d^)2 

•> "Ar^^djdj^ - ^(Cjdj + e^d^^ + 1) + 

- i[c^(dj + dp - 2(1 + Kpj)(l + fjdpi 


2T 


fe(dj 


‘^k> 


(38) 


wheie T is given by Equation (28). Since ^ shows 

that Equation (30) is redundant if Equation (31) is enforced. Consequently, 
the satisfaction of Equations (33) and (35) constitutes the necessary and 
sufficient condition for unconditional stability. 

Expliciting Equation (35) requires some further algebraic treatment 
which is presented in Appendices C and D, From this, if one lets 


w ■ (0 - ]/2)/9’, the following stability conditions result: 

0(' e " *^6'^ - ^i ®(‘e (39a) 

~2 (t'e - 2u) < (39b) 

for the case of second-order smoothing (see Appendix C), and 

20(2ej, - t^ucg) i i 2e4iT^ (40a) 

2e(e-t)iS <^0b) 


for the case of fourth-order smoothing (see Appendix D). The corresponding 
donuilns of 'hmc»>nditional stability** are shown on Figure 1, Note that in 



(a) Second-order smoothing applied explicitly. 


e ■ ISe, 



(b) Fourth-order smoothing applied explicitly. 

Figure 1.- Domain of unconditional stability for periodic boundary 

conditions. 
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case fourth-order 8K'‘««-»-tng Is applied, this domain Is bounded (Figure lb). 
For the trapezoldal-tlme-dlfferenclng scheme, the domain reduces to the 
line * Cg/2 In case second-ordor smoothing Is applied, and collapses 
to the origin In case fourth-order smoothing is applied. These conservative 
results conflict with the authors' numerical experience (see Chapter III). 
This was attributed to the assumption of periodic boundary conditions which 
results In an Inadequate model relaxation problem, as mentioned at the begin- 
ning of this section. For this reason, the case of specified boundary data 
Is examined In the next section. 

3. The case of specified boundary data 

When the solution u** Is specified at the boundaries, the forward 
shift operator, Trld (0,0,1), and the backward shift operator, Trld (1,0,0), 
are no longer Inverse of one another, and cannot be simultaneously diagon- 
alized as they could for periodic boundary conditions. (In fact, they are 
both singular and In Jordan canonical form, or the transpose of It, with 
all the eigenvalues equal to zero.) As a consequence. In Equation (11), 
matrices with the same subscript (x or y) cannot be expressed as linear 
combinations of (positive or negative) Integer powers of a unique (shift) 
operator, and cannot In general be simultaneously diagonalized. This Is 
true. In general, for the matrices and together, and Ay and By 
together. (If second-order nnoothlng Is applied, a case of exception occurs 
when - 0c^.) However, approximate commutation of the right- and left- 
hand sides of Equation (11) occurs when the coefficients of the smoothing 
terms are either very small or very large compared to the coefficients of 

the convective derivative operators and Cy. This situation corresponds 
to 



I ♦ f ♦ I 


— <^ 1 or >> 1 


where c is either or and v is either or Vy. 

In practice, implicit smoothing Is Introduced In an attempt to keep 
the coefficient directly proportional to At, and still maintain uncon- 

ditional i<tabillty. If one assumes a priori that this is possible, and 


Ce ' ag At 


Ci ■ 0^6 At 


for some constants a,, and Equation (41) becomes: 

Av AV 

— and ^ « 1 or >> 1 (43) 

The mesh spacing parameters Ax and Ay can certainly be considered as very 
small and so are Ax/a and Ay/b, In general. However, if this inod:'l prob- 
lem Ic of any relevance for the Euler equations, the wave speeds a and b 
should play the roles of the eigenvalues of the Jacobian matrices A and B. 
Scnw’ of these eigenvalues can eventually become very small in some regions 
of a tranfionic flow field (see Chapter IV), so that both limits in Equa- 
tion (43) are of interest. If one makes the assumpjlon that these extreme 
Situations ((1) Ax/a and Ay/b very small, and (2) Ax/a and Ay/b very 
lar^.o) ar<' those that produce the binding conditions for stability, one is 
tempted to analyze the as>nnptotlc properties of the algorithm In these two 
limits. This Is precisely what Is done in the rexnalning part of this 
sect Ion. 


^To bring a theoretical support to this assumption Is precisely the 
motivation for this analysis. 


L 
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ronBider first the case where Ax/a and Ay/b are both very small. 

This case will be referred to as the case of large Courant numbers and 

Vy). Choose X and Y to be the tranr format ions that diagonalize the 

matrices C„ and C„ respectively. These are given bv (see Appendix B) : 

A y 


X^j - .^2/(J + 1) 1® sin mOj 


/2/(K 4 1) 1™ sin me 


‘mk 


k ; 


(44) 


where now t) ■ Jn/(J 4 1) (J •» 1,2, . . .,J), and 0 « kr/(K 4 1) 

J K 

(k » 1,2, ♦ . .»K). As for the case of periodic boundary conditions, these 
transformations are unitary (C and C are skew-symmetrlr ) , so that: 

XV 


x;;,] - = /2/(j 4 D(-i)-’ sin ,ie.^ 

y;j* - Y*,, - »'2?(K 4 l)(-i)‘‘ sin ke 


‘mk 


m ; 


(45) 


As a result of this choice, the matrices C and can bo written a- 

A y 


follows: 


wheie - Dlag(cj) and Ky 
C|^ • cos On their part, 
become: 


- X(1K^)X*’ ' 

> 

Cy - Y(lKy)Y-' 

«• Dlag(cj^), In which now, c^ - cos 6j 
the matrices and in Equation 


(46) 


and 

(18) 


and the matrices Ay and By are given by similar equations. In these equa- 
tions, matricea proportional to or Vy are now considered as principal 
parts, and the oth'r matrices as perturbations. Recall that as a particular 
case of application of a general result of perturbation theory (see, e.p., 




> 1 ' 




(17)), the first-order perl urbnt ions on the elgenvnluv->s of a dlaKO'>al 
matrix (with distinct eigenvalues) are simply the diagonal “lements of the 
matrix by which It Is perturbed. In view ; Equations (17) and (i*?). It 
appears that the off-dlagonal elements of A are themselves flrst-ord'ir 
perturaatlons and thus contribute to the eigenvalues of A by terms that 
are at least se-ond-order perturbations. Su.h perturbations are neglected 
In the remaining part of this derivation. In this appi *matlon, the 
matrices ami 6^^. for example, become diagonal matrices with eigenvalues 


glv« n by 


aj - 1 + dv^lcj + t 
bj - v^lcj + r^dj 


where dj and dj are defined by 


^.1 ■ 


<1; - 


and terT.8 of order (c/v^)- are neglected with respect to one. But Equa- 
tion (48) is analogous to Equation (22). Thus, unconditional stability Is 
enforced bv Equ.it Ion (IS), provided dj , d^. d ' . and d' are replaced bv 
*^J* *^k’ rp'^ppctlvely, which act as ’V^fectlve eigenvalues" of 

the smoothing operators at large Courant numbers. These effective eigen- 
values are evaluated In Appendix E. In particular 

\ " 2 (SO) 

(which is the average value of d^ or d^); making the corresponding substi- 
tutions In Equation (3S) yields, ./ter some simplifications: 


*p«'-'(dj + dj^) - 4P(2r^ + 1) + 2 V 0 
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or 


Y0 


- 7 ) 


(51) 


where again u * (6 - l/2)/9^, and y - (l/4)Max(dI d'). In particular, 


if second-order smoothing is applied, d^ » d^ 
and Equation (51) becomes: 




d,^ = d^ = 2, so that > = 1 


(52) 


If instead, fourth-order smoothing is applied, Max (dp and Max(d|) converge 
to 8 in the limit of a mesh refinement (see Appendix E)„ so that y 4 and 
Eqtuitlon ( 51) becomes 

'i - ('e - !) <”) 

One can observe an analogy between Equations (52) and (53) and Equations (39b) 
and (40b). 

Consider now the reverse situation where the Courant numbers v and v 

X y 

(alLern&tt i.y the wave speeds a and b) are small compared to the coeffi- 
cients and of the smoothing teims. For this case, the matrices 
X and Y appearing in Equation (17), are chosen to be the orthogonal 
matrices C and n that diagonalize the (real symmetric) smoothing operators 
Djj and Dy (and also 0^ and Dp, respectively. (The matrix C is explicited 
in Appendix B. Tt Is found symmetric, but this property is not used here.) 

In this way, the smoothing operators D and D (and D' and D'), which pro- 
ducf the principal part of A, are represented by diagonal matrices, while 
the convective derivative operators and VyCy, now sought as perturba- 

tions, arc represented by the following similar matrices: 

-t. 


^x^x 


V Cx' 


* ^’yH CyH 


( 54 ) 


' I 





m 






:s 


t 


It stH'oiui-’orJor |u»rt urb^tt Ions ow iho oigeiuMluos of \ aro noglorteti, only 
tho diagonal olt'nuMUs of and uood to bo rofalned, Thos«' are /.oio 
boi\uiMo t j. and C’y , and oonsequont I v and (*^, are skow-svironot r iv' • llonco, 

lirst-ordor analysis and /t'rot b-ordor analysis of this oaso produoo t ho 
sank rt\sul t . Tho lattor ono ot»nslsts ot treating the oaso o»f pure diffusion 
for whloli Equation (2 0 applies If, In tho definitions tliat follow this 
oquat Ion (Kquatlons (24) and (2S)), tho parameters, Cj and o, are sot equal 
to :oro, and If 0 mid b are dot Inod as In Fquatlon (4*^). Since those 
are tho only nKvli f Icat ions to l>rtng to the analysis developed for periodic 
boundary conditions, tho stability rendition is given bv Kquation (30), or 
oqui valent \ v : 


2(1 -i 


I'V 


.1 + 


■1‘V 



(Vi) 


Enforcing that Kquation (SS) bo satisfiod for all values of 1 and k 
resulted in the following condition: 



tor tlu' Ciiso whore second-order smootlilng Is applied, and In 


(Sb) 



(37) 


tor the case whore fourth-order smoothing is applied, (The derivation of 
Ki|uations (SeO uid (S/) Is given In Appt'ndix K). 


In ci'MchisIon, an approxliiuate definition of the domiiln of uncond 1 1 lonal 


staMIltv stunilvl be c>btained bv combining t lu' condition!^ given In Kqua- 
tions (S2N and (Sb) for the case where second-order smocUhlng is applied, 
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and in Equa^.lons (53) and (57) for the case where fourth-order smoothing is 
applied* These domains i^re shovm in Figure 2 tor the trapezoidal time dif- 
ferencing method (9 « 1/2), and in Figure 3 for the Euler implicit method 
(0 1), On the latter figure, the domains that have been previously 

obtained assuming periodic boundary conditions, are reproduced for compari- 
son. It appeals that in assuming specified boundary data instead of 
periodicity, results in less stringent stability limitations. 

The key result of this analysis is that it suggests that the domain of 
unconditional stability in the (e^, e^)-plane is unbounded, allowing the use 
of arbitrary values of At and provided is maintained sufficiently 

large. (This will be demonstrated in the next chapter by various numerical 
experiments that were conducted on a more complex problem.) In practice, it 
should be sufficient to let 


ece 2 


or 2 


(58) 


when elthei second-order or fourtVi-order smoothing is employed. If now, 

is kept directly proportional to At, the inconsistency of the algorithm 
is removed, and this, theorelically , without violation of the property of 
unconditional stability, provided Equation (58) is enforced. This was not 
possible with the original formulation (Equation (2)). 


C. Application: Sequence of Parameters 

lu Section IIB above, it was shown that in the modified differencing 
scheme, Equation (3), the numerical dissipation could be kept directly pro- 
portional to At. If this is done, the steady-state solution Is independent 
of At, and much larger values of At can be taken without triggering 
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(a) Second-order smoothing applied explicitly. 

(b) Fourth-order smoothing applied explicitly. 

Figure 2.- Domain of unconditional stability of the trapezoidal time 

differencing scheme. 
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(a) Second-order smoothing applied explicitly. 
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(b) Fourth-order smoothing applied explicitly. 

3.- Domain of unconditional stability of the Fuler implicit method. 
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1 . 
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nonlinear Instability. Large values of At, as well as a sequence of a 
small CO a largo At might thus be used to accelerate steady-state 
convergence. In this section, some motivations for using this technique 
are given. Consider again the simple model two-dimensional first-order 
wave equation. Ivquatlon (10). If the solution u is specified at the 
boundaries, and if numerical dissipation is not applied, the eigenvalues of 
the iteration nuitrlx L can be obtained from an equation similar to 
Equation (23), and are given by: 

1 - Cl 6j cos 0^^ + 1(0 - l)('>j{ cos 6j + Vy cos 0y^) 

^jk * 1 - v,-v„ cos 0. cos ", + ’ (v cos 0. + v cos 6.) 

xyj .^jy^ 

where * ah/Ax and ■ bh/Ay. The h^-term which appears in the real 

parts of both numerator and denominator of is the cross term that 

results from the approximate factorization of the left-hand side of the 

difference equation. 

It is directly apparent that for the trapezoidal time differencing 
method (0 * 1/2), the modulus of is exactly equal to 1, whether 

approximate factorization is used or not. This means that no dissipative 
mechanism exists to permit the steady-state convergence of this method, 
unless smoothlni' is applied or boundary conditions are modified. 

Consider now, the Euler Implicit method (0 • 1), and rewrite X.. as 

3 ^ 

follows : 


X 


jk 


1 + 1<|. 


jk 


(60) 


where 


h((a/Ax)cos 0j + (b/Ay)cos 6j^] 
Jk 1 - h2 (a/Ax) (b/Ay)cos 0. cos 0j^ 


( 61 ) 


r 


I 


i 




.he.- 


■i^n A .Afc. 
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Without approxlinate factorization, the h^-tenn in the expression of (p 

J ^ 

would disappear. Then would be proportional to h, and 1^,. | 

would decrease with increasing h. In this case, the larger the time-step, 
the faster the ct-ivergence would be. However^ a different conclusion can be 
drawn If approximate factorization Is used. For chat case, for given wave 
speeds (a and b), and frequency parameters (0 and 0 ), the modulus of \ . 

j ^ j ^ 

achieves a minimum when Ujj^l maximum, and this occurs, when 


h » 


Ax 


Ay 


a cos 0, 


b cos 0, 


1/2 


(62) 


Thlt shows that there exists an optimum time-step parameter h which not 
only depends on the wave speeds, a and b, and the mesh spacing parameters, 

Ax and Ay, but ilso oi. the frequency parameters, 0 and 0 . For small 

j 

values of cos 'j and cos (Interpreted as low frequencies), or for small 
wave speeds a .ind b, a large time-step is desirable, as anticipated. How- 
evei , for values of |cos 0^ | and |cos 0^^ [ of order 1 (Interpreted as high 
frequencies), the Coinant numbers ah/Ax and bh/Ay should themselves be ol 
order I for a ripid reduction of the residuals, To Illustrate this, the 
range of values of cos 0 and cos 0 for which |^,, | i C (a con.?tant 
taken to be 0.9>), Is represented in Figure 4 assuming Ax/n ■ .* . for 

different values of the Courant number v ■ ah/Ax - bh/Av, On this figure, the 
corners (|cos 0|] « 1 or |cos ('. | ■ 1) are Interpreted as high frequency 
regions, and the neighborhood of the lines cos 6, ■ 0 and cos 0. ■ 0 as 
low frequency ri*glons. The circles represent the values actually achieved 
by cos Gj and cos for a mesh containing 15x15 Interior grid points. 

The domain ^ consists of two strips that converge tow.ards the low 

frequency region when At Increases. 
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From this, one concludes that a large time-step should be efficient at 
low frequencies, but also, that a sequence from a small to a large time-step 
should be efficient by nc*. privileging any particular frequency band. 

These concepts have served as a guideline for the numerical experimen- 
tation of the next chapter. 


1 
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TH. NUMERICAL EXPERIMENTATION ON STEADY- STATE CONVERGENCE 

In this ch.ipfer. a model problem governed by the Euler equations is 
solved numerically to compare the iterative convergence properties of the 
modified algorithm, given by Equation (3), to those of the base algorithm, 
given by Equation (2). 


A. Model Problem 

A model transonic flow problem was selected to test the convergence of 
the modified differencing. In the past the transonic flow about a nonlift- 
ing biconvex airfoil with linearized boundary conditions has served as the 
prototype problem for relaxation algorithms and so this problem was used 
hero. A variable grid with clustering was used to resolve flow-field 
gradients (see Figure 5), but the equations are solved on a uniform trans- 
form plane by introducing simple stretching transforms. 

The solution procedure is as follows. The values of the conservative 
variables at interior points are first advanced from some starting solution, 
using either Equation (2) or Equation (3) with h « At. (The Euler implicit 
roetliod (0-1) Is preferred here, to the trapezoidal time differencing 
method (0 - 1/2) which is nondlsslpatlve (see Section IIC), because the 
emphasis, in these numerical tests is on steady-state efficiency.) Then, 
verv simple boundary tonditlons are applied. Free-stream conditions are 
enforced at the inflow and upper boundaries. Along the body, the y compo- 
nent of vt'locity is obtained from thin airfoil theory: 

V - lU(dy/dx)^ (63) 

where Is the fre<»-stream velocity, and (dy/dx)g the body slope which 

Is a specified function of x. All other unknowns are obtained by 
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(zeroth-order) i*xtrapolatlons. Higher-order boundary conditions improve 
accuracy, but are deliberately avoided in this study because they can sig- 
nificantly degrade the stability and convergence properties that one would 
prefer to isolate. 


B. Results 

Results for a lO-percent-thick biconvex airfoil at M<„ ■ 0.84 are 
shown in Figure 6, and compared to a potential solution by Holst [8). It 
should be noted that a coarse grid and simplified boundary conditions have 
been used in order to test a variety of parameters. Much better solution 
accuracy is obtained by grid refinement and use of more accurate boundary 
conditions. Detailed solutions of this nature are available in [4). 

The solution shown In Figure t> was obtained using cither Fquatlon (2) 
with a nondimenslonal At ■ 0.03 and » 0.03, or Equation (3) with a 
nondimensional At • 0.38, 0.38, and These values wore 

each found to be close to optimum l>y a trial and error process. The con- 
vergence histories for both cases are shown in Figure 7 where root-mean- 
square residual error as well as the average difference between the con- 
verp.ed and intermediate Cp distributions are indicated. Recall that the 
boundary conditions are applied in an expllcit-like manner, which is 
expected to slow the more raplilly converging case, that is. Equation (3) 
mon significantly th.nn Equation (.!), which uses more time-steps. Figure 7 
shows that the modified differencing converges to steady state about 8 times 
faster than the original scheme. This experiment tends to verify the con- 
clunlon diawn from the model problem - a large value of At c.in be effec- 
tive in ai hleving more rapid steadv-state convergence. 
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Figure <1.- 


Converged pressure distribution along the airfoil. 
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c 7.- Effect of using numerical dissipation implicitly (Euler implicit 
diffeiencing scheme wi*-h = At). 
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In these experiments, the starting solution was taken tJ be the one 
obtained after 25 applications of the base algorithm (At ” ■= 0.0., 

« 0). with all the properties Initialized to their free-.‘?tream values 
and gradually Introducing the body by increasing with time the body slope 
(dy/dx)g from 0 to its correct value. In this way, impulsive starts were 
avoided. 

It must be noted that the ratio of to can significantly influ- 
ence the conver;^ence rate. It uas verified that for = At, and a single 

optimized time-step, this ratio could optimally be set equal to 2 (for the 
Euler implicit method), as Indicated in Figure 8. For larger values of 
this ratio, the dissipation term added implicitly excessively stabilizes the 
transient behavior of the solution. For smaller values of this ratio, the 
coelficlent tg, and consequently At = Cg, must be reduced for stability 
(see Figure 3b), and this reduces the rate of convergence. 

The effect of using a sequence of At is Indicated in Figure 9. In 
order to simplify the optimization of the sequence, the fol owing formula 
was used 

^tn - At, + (S^) - •'":) 

where n 1,2,. . .,N for a vvcle of N time-steps, e * 2 in most experi- 
ments, ami At^ and Atj^ were optimized. The data show that a sequence of 
At Is efiective but not as lauch as one might expect. The sequence of 6 .^t 
Is about 10 times more effective in steady-state convergence than the 
orij;lnal scheme. The data shown are for optimum values of i p and At, 

In (omparison to the same scheme based on using a single optimized time-step, 
a 3(*quenc4 of parameti*rs saves 50 to 75 time-steps out of 150 to 250 time-steps. 
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Figure 8.- P^ftect of the ratio of to ' e (Euler Implicit differencing 

scheme with = At), 





SUIJ/ 



ITERATIONS, n 


Figure 9.“ Effect of tislng a sequence of time-steps (Euler Implicit differ 
enclng scheme wl*'h e_ - At and Cj * 2tg). 
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depending on what solution tolerance is desired. In this case, plottable 
accuracy Is actually achieved after 150 time-steps with the most effective 
sequence. Again the model problem preaicts the iterative convergence proper- 
ties of the more complicated flow, and the use of a sequence of time-steps is 
also an effective way to accelerate steady-state convergence. In these tests, 
the optiraum values of At^ and At^ were found to correspond to a limit of 
stability. It was also noted, that at this limit, the average value of At 
for a cycle of N time-steps is the same for all the cases shown in 

Figure 9. 

Note that more sophisticated procedures for controlling various param- 
eters would lead to better convergence rates. Tn particular, it was 
observed that a more rapid convergence could be obtained (for this problem) 
by setting =* 1.02z^ (instead of * 2e^), » At and choosing a 

sequence of time-steps that Includes one or two that are sufficiently large 
for (c^, i to fall in the unstable range. It also appears that the 
operational range fAt^, Atj^] should be optimized with the solution itself, 
tha! is, with the iteration counter n. It is most likely that these would 
produce better Lmprovt^ments. Nevertheless, they have been avoided here 
because of their lack of simplicity and generality. 

Sensitivity in reta of convergence to nonop tiinality is weaker if N 
is large. For example, Figure 10 shows that for a cycle of 6 time-steps, 
if Atj and Atg are set equal to half of their optimum values. At* and At*, 
the algorithm, over the first 300 steps, loses less than 20 percent in rate 
of (onverjence, and remains as efficient as It is fo'- a single optimized 
time-step. 
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The tlfect of varying the exponent e for fixed At^ and At^- Is 
Indicated n Figure 11* As e decreases, the awrage Lime-step increases 
and so does the rate of convergence until nonlinear instability occurs 
(e 1). 

Sensitivity in «*ate of convergence to free-stream Mach number was 

also studied as indicated in Figure 12. Three cases were computed using 
the same sequence of time-steps; that is, the sequence was not optimized 
for Mp,,. The data indicate that the implementation of implicit smoothing, 
and the use of large time-steps extend to subsonic and supersonic regimes 
as well as transonic regime. 

The influence of the boundary conditions on the rate of convergence 
was also investigated. In this test, a sequence of six time-steps given by 
At^ ■ O.OS, 0.2, 0.45, 0.8, 1.25, and 1.8 was used, and « At « e^/2. 
This sequence was not optimal, but this is not believed to have had any 
importance. At first, a fully converged solution was obtained. The start- 
ing solution was then constructed by increasing by 5 percent the converged 
solution at interior points. The rate at which this disturbance could be 
eliminated, for some given boundary conditions, was then evaluated by com- 
puting the following estimate tor tV»e spectral radius of the Iteration 
matrix: 

where RMS^ is the root:-mean-square residual error (the right-hand side of 
Equiition (3)) after n applications of the algorithm. An estimate for the 
numlier of time-steps required for a reduction of the residual errors 
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by i\ factor of 10 was then computed according to the formula: 
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Figure 11.- Rate of convergence for various se .^..ces of time-steps (Euler 
Implicit differencing scheme with ■ At and tj ■ 2e ). 
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Flj;ure Effect of the free stream Mach number (Euler Implicit differ- 

encing, scheme with ■» At and Ej » 2e ). 
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n^Q “ l/colog2o(p) 

Four tests were made. In these, the inflow and upper boundary values were 
fixed to free-ntream conditions as normally done. However, the body bound- 
ary and outflow boundary values were either fixed to their converged values 
or variable, as permitted by the regular extrapolation procedure. The 
values of p and n^Q obtained in the four cases are collected in Table 1. 
The results show that boundary conditions have a very strong effect on con- 
vcrj^ence properties. By fixing the body-boundary values (although this is 
impractical since it requires prior knowlet^ge of the solution), the rate 
convergence doubled from what it was in the regular procedure. This favor- 
able effect is even stronger if Instead the outflow boundary values are 
fixed. This case is more practical for transonic flow applications where 
the properties at the outflow boundary can be fixed to free-stream values 
without significant degradation of the solution accuracy. If now all four 
boundaries are fixed, an improvement in rate of convergence by a factor of 
4 is observed. This experiment indicates the strong dependence of the 
Iterative convergence properties of the algorithm on boundary conditions, 
and opens a possible area of investigation for future work. 

Table 1. Influence of the boundary conditions on the rate of convergence 


Outflow boundary 
values 

Body boundary 
values 

P 

"lO 

Extrapolated^ 

Extrapolated® 

0. J8800 

191 

Extrapolated 

Fixed 

0.97633 

96 

Fixed 

Extrapolated 

0.96881 

73 

Fixed 

Fixed 

0.95064 

45 


Regular procedure. 
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Finally, the effect of using only second-order smoothing, explicitly 
as well as Implicitly, was Invest iy;ated. For this test was set equa** 

to At/2. The base algorithm = 0, and « At/2) was found to operate 
optimally with At * 0.22. If insiead, a sequence of time-steps is 
employed an improvement in rati* of convergence by a factor of 2 or so is 
achieved, as shown in Figure I'K This test also shows that the use of 
second-order smoothing considerably Increases the rate of convergence of the 
regular algorithm (c^ * 0) itself. However, unacceptable losses In accuracy 
occur if this type of artificial dissipation is employed to calculate a 
flow field with a large change of gradient in the solution. Even in the 
sim]>le biconvex airfoil calculation considered, the solution at the leading 
and trailing edges is noticeably degraded, although the shock wave is still 
adequately resolved. 

One concludes in general that large At is very effective and that use 
of a sequence of At can be perhaps twice as good. The algorithm is not 
overly sensitive to nonoptimum features. However, better rates of improve- 
ment seem possible (e.g., the added effectiveness when only second-order 
dis^;lpation is used, better boundary condition procedures). 

One remarks that some ideas have proved effective as well in more com- 
plex flow calculations [4-6]. However, the additional sensitivity of the 
mon* complex flows to nonlinear Instability forces the use of much smaller 
t inn-steps. Consequently, the improvement in rate of convergence is much 
les.s “ typically a factor of 3 or 4 over the base algorithm. 
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IV. ON THE EFFECT OF SPACIAL VARIATION OF THE JACOBIAN MATRICES 


In the cossputatic-a of iransonic as well as subsonic flows, the eigen- 
values of the Jacobian matrices A ■ [SE/ctq] and B « (3F/3qj are real, 
and of mixed sign. This is the reason for adopting a time-dependent 
approach, combined with the use of central space differencing which, in 
principle, produces purely imaginary eigenvalues for the convective 
derivative operators 

Cjj « (2Ax)6j^A 

(64) 

Cy =• (2Ay)cS^B 

The unconditional stability of the implicit algorithm, derived in Chapter II 
for a scalar, linear model equation, relies crucially on this properly. In 
this chapter, the possibility of breakdown of this property for the Euler 
equations due to variable coefficients is examined. 


A. Analysis 

In this analysis, Che matrix Is considered, in particular. A mesh 

containing, J>^K interior grid points is assumed, so chat the dimension of 
the matrix is (J^K'<p)^ for p dependent variables (p *= 4 for twj~ 

dimt*nslonal flows), lor convenience, it is here assumed that the 
components of the solution vector q are ordered as follows: 

/tt ttt V tt tvC 

^ ^21* ' * *’ ^ ^J'’’ * * 


whe}*e q contains (he p dependent variables evalratod at (x., y, ) . 

J ‘ k 

Then, making the simplifying assumption that Ihe solution vector q 
Is jipecifiied at the bcaindarles, permits us to write the matrix as follows: 
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- BDiag(C^ ,^^) (k - 1,2,. . . ,K) (66) 

where C„ i, is the following Jp'^Jp matrix: 

tv 

<^x.k - “• ''j + l,k> <J ■ = ■'> 

In vhich Is the p^p Jacobian matrix A evaluated at (Xp v^^). 

deeply, if A vas some symrat^trlc constant matrix, the matrix 
would be real skew- symmetric and would trdeed have purelv Imaginary eigen- 
v'alues. However, one may question whether this property carries over to 
the general case where A Is *1 nonsymmetrlc matrix subject to appre- 

ciable variations from point to point, due to the nonuniformity of the mesh 
as well as the solution itself. The purpose of this analysis is precisely 
to bring some information about this question. 

It first appears from Equations (6b") and (67) that the eigenvalues of 
the matrix are obtained by collecting those of the matrices ^ 

togc'ther. For this r<*ason, only one of these matrices will now be consid- 
ered, with the subscript k omitted in what follows. 

It is known that tlie Jacobian matrix A tor the Euler equations can 
be tl iagonalized by a real transformation T, so that: 


A »= TAT'^ (68) 

whei e A *= Diag(a”^) and m » 1,2,. . .,p. Explicit expressions for T and 
A I an be found in [IH]. For example, for a two-dimensional flow, if 
Cartesian coordinates are used: * a-' » u, ~ u + c, and » u - c 

where u is the x component of the velocity vector and c is the local 
spet'd of sound. For the same case, the eigenvalues of the Jacobian matrix 
B are b * b* * v , ^ - v -f c , and b"* « v - c, whore v is the v com- 
pont'nt of the velocity vector. 


t 


\ 

\ 


k 





I t 


I 


I 
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Now, construct the following matrix 

BDiagd^T^) (69) 

and perform on the matrix Cj^(truly some k) the following sim- 

ilarity transformation: 

- BDlagl(-l)^TT']BTtI<J(-A^_j, 0, )BDlag(i^Tj ) (70) 

* ia 


where 


a . BTrid(IJlTj_,A^.,,O.IJ>Tj^,Aj^,) 


(71) 


It is desirable that all the eigenvalues of the matrix o be real and for 
all those of the matrix to be purely imaginary. 

Note that If the flow variables are continuous, the matrices i 

and 

j j*+ 1 depart from the identity matrix only bv terms of 0(Ax). For 
this reason, one expects the eigenvalues of the matrix o to be well repre- 
sented by those of the fc Hewing matrix: 

» BTrid(Aj_j, 0, - e (72) 


In Tuaking this approximation, one assumes the effect of variable coeffi- 
cients to consist vripiarlly ot ^he variat'^^^n of the eigenvalues of tliO 
matrix A^ rather than the variation of its eigenvectors. This assumption 
is made ht're, and the next step consists of rearranging the rows and tne 
columns oi the matrix a’ to collect the eigenvalues a.^ (for given m) 
tog(‘ther. More precisely, for some nonsingular matrix P, which corresponds 
to a product of permutations, the matrix 




*r 


(73) 
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which is s-'nilar to the matrix o', becomes 

a" = BDiag(r*") (m»l,2,...,p) (’A) 


where 


Trld(a”^^, 0, 


'J = 1,2,. . .,J) 


(75) 


Hence, the eigenvalues of the matrices o’ and o" are obtained by collect- 
ing those v f the matrices r”'. These eigenvalues must be real for those of 
the matrix to be purely imaginary. In this way, the analysis is 

reduced to the one of p independent scalar problems. Thus, one is lead 
to exaMne the eigenvalues of the matrix for a particular value of m, 

and to omit this superscript in what follows. This is done in Appendix G, 
whose main results arv repeated here without derivation. 

It turns out, that for all the eigenvalues of the matrix F to be 
real, it suffices that the following condition holds: 

° * -tJ - 1) (76) 

This condition is met either when each eigenvalue aj of the Jacobian matrix 
A lias the same sign at all the grid points, or when it does change sign at 
one or more grid points but vanishes exactly at one or more grid points 
before changing sign. The condition given in Equation (76) is not, however, 
necessary for all the eigenvalues of F to be real. Nevertheless, this 
favorable result is most unlikely to be true If this condition is violated. 

To see this, another result of Appendix G is recalled. For this, define 
sequences of coef ficlt'nts ^rid by the following recurrence 

formulas: 


(v) „(v-l) ^ (v) 

t * 6 -Fa .a , a 

nrfl m ?m+2 m 

i^; - H- a ^ a ^ 

mfl irrfl ?tn+2 ?ir&3 


cm 
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wheie m is a natural Integor and v « 0,1,2,. . . ,m, and the followin;; 
conventions are adopted: 


. 1 

m m 

- 0 

m 


(78) 


These definitions being made, it turns out that a necessary condition for 
ail J eigenvalues of the mati T to be real is that the coefllclcnt 

(v) (v) 

a , In case J • 2m, or 8 , in cose ^ ■ 2m + 1, be positive for 

m m 

V * 0,1,2,. . . ,m. It is easy to calculate in particular and 

in m 

which ar'' given by: 

. (ol 


m 


a 


• a. 


m 


. (o) 


a a . . .a I — + — + . , , + ] 

' ■ V’l 


(79) 


In genera], and f ^ are polynomials of degree 2(m - v) of the coef- 

ficients i\y and It is most unlikely that they all will remain positive if 
the condit Ion given in Equation (7b) is violated at one or more grid points. 

In particular, situations where alternatelv Is negative should 

m * m 

be comnx^n. If this happens, the nvitrlx V has an odd number of pairs of 
purelv Imaglimrv eigenvalues, sav i r^, and -i (C - 1,2,. . .,2g + 1). 

To these correspond the real eigenvalues -r , and r, for the matrix (' , 
hall of which are negaMve. For trapezoidal time differencing end arbitrary 
tlmt-slep, or Euler implicit differencing and sufticiently large time-step, 
tiu'se real negai ive eigenvalues produce numerical instability unless thev 
are I'alanced lo' a suttlclent pivsltlvo contribution coming from the smoothing 
opoiators. The;o unfa'^orable I'igenvalues should, however, be of snuill 
modul I If oiU' assumes Muit they are essentially determined by the entries 
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of the matrix in the neighborhood of the point where the alternation 

of sign occurs (weak destabilizing effect). However, since the matrix 
appears in the difference equation multiplied by At, the coefficients of 
the smoothing operators should themselves be kept proportional to At as 
required for consistency. 

In conclusion, it appears that a particular form of Instability due to 
variable coefficients may be triggered If central space differencing is 
used at a point where one of the eigenvalues of the Jacobian matrix A or B 
changes sign. However, use of numerical dissipation proport lonal]y to At 
and in sufficient amount, should remedy this type of instability. 

In the following section, some numerical examples of this phenomenon 
are presented for some scalar model problems. 


B. Numerical Experiments on Scalar Model 
Equations with Variable Coefficients 
In the numerical experiments, the trapezoidal time differencing method 
was used because this method is neutrally stable, that is nondissipatlve, 
for scalar linear problems with constant coefficients. In this way stabil- 
ity problems due to variable coefficients could be isolated more easily. In 
all tne cases, scalar functions of the only two Independent variables x and 
t were considered. 

In the first test , the following problem was solved 

Uj. 4- [a(x)(u - Dlx - - 0 (-1 < x < 1) \ 


(80) 


where the wave speed 


u(x,o) - Uq(x) [• 
a(x), chosen to be 


2 exp(-5x^)] J 


a(x) - 4x/(l + 27x**) 



» I I 
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had the si^n of x. For r * 0, this probli’m Is a rarefaction wave problem. 
The characteristic curve In the fx,r) plane has the slope dt/d>: = l/. (x) 
and is pointing outward cf the donviin of integration at the endpoints 
X = ♦!. For this reason, specifying the solution at these boundary points 
would here be improper. Instead, in the numerical computation, (zeroth- 
order) extrapolation was used at these points. In this way, some snvill 
positive terms were introduced in the diagonal of the matrix at the 

upper-left and lower-right corners. The numerical computation of the' 
eigenvaiuis of the matrix confirmed that these terms produce some 

positive contributions to the real parts of the eigenvalues (compared to 
the case where the solution is specified at the boundaries), and thus have 
a favorable stabilizing effect (outflow of residual errors). Despite' this, 
with c » 0, the trapezoidal time differencing method was found unstable 
when using a mesh with grid points located at Xj >= (J - 16.5)/15 
(j ' 1,2,. . .,32), so that a (x j^^)a (x^ ^) ^ 0. This instability renujlned 
for value.', of les.s than 6 '< 10“**, or so. For larger values of r, the 
numerical solution remained bounded and in fact convergent, at all the grid 
points, ti' the exact .•^teady-st.ite solution of this pr obi em which is 
u(x,'") = ,1. This steady-state solution was not altered by the smootlilng 
terms in this ideal case where u(x,"') is constant in x. It was also 
verified that a stable but not convergent (to steady-state) algorithm was 
obtained lor e « 0 by locating the grid pi Ints at x^ =• (x - 16)/IS 
(J ' 1.2,. . .,31) so that positive values of a (xj ) were separated from 

negative i>nes by a true zero. 
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1 


( 81 ) 


In the raoalnlng tests, the following class of problems was considered; 

Uj + [<|'(u,x)]^ - eujjx “ [4>(Uoo(x).x)l^ (0 < x < 1) 

u(0,t) = u<^(0) 

u(l,t) * u„(l) 

u(x,0) =• Uq(x) (specified) 

in which Uo„(x) - exp(-5x^), and the functional forms of 4> and varied 
from case to case. The particular form of this equation was chosen antici- 
pating that for sufficiently small c, the stationary solution of this 
problem would approximate the lunction u„(x). This function was chosen 
rather arbitrarily but nonuniform so that ^ 0 at the steady state. 

Mor»'over, in the numerical computation, the term appearing on the right- 
hand side of the differential equation, that is the source term, was cen- 
trally dll ferenced in the same way that the corresponding term of the left- 
hand side of the equa Ion, that is, the flux term. In this way, the 
sn»othing term was entirely responsible for the discrepancies 

between u(x,«>) and Uj..(x). 

Although, to the author's knowledge, specifying the solutions at the 
boundaries always leads to a well-posed problem for c > 0 (assuming smooth 
data), this is not necessarily the case for ■ 0, as mentioned prevlouslv. 
In particular, for the latter case, the following inequalities 


3u 


ii 

Du 


> 0 


x«o 


< 0 


( 82 ) 

( 83 ) 


x*l 
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should hold for the characteristic curve to point inward the domain at the 
boundary points. In all the rases that follow Equation (82) applied, but 
not necessarily Equation (83). This question will be discussed for specific 
examples. 

Ten experiments were conducted on linear test equations that were 
obtained by letting the flux function be of the form 4>(u,x) *= a(x)u. 
These experiments are defined in Table 2. For the first three cases, a(x) 
was chosen strictly positive. For this reason, the implicit algorithm was 
found stable. Kowever, adding artificial dissipation was found necessary 
to obtain a steady-state solution. For a rather small value of t. (Test 
No. 2) the numerical solution is very accurate as shown in Figure 14. For 
this case, the slightly improper boundary condition u(l,t) *= it^,(l), does 
not disrupt the stability of the algorithm (even for c « 0), and does not 
seem to degrade the solution accuracy significantly (for c ^ 0). This 


Table 

2. Numeric 

:'al experiments for 

linear 

test equations 

Test No, 

a(x) 

Uo<x) 

L' 

J 

Comments 

1 

3.5-Hx 

exp (-5x) 

0 

52 

Neutral iv stable 

o 

3.5-3X 

exp (“5x) 

0.025 

52 

Convergent 

3 

3.5-3X 

exp(-Sx) 

0.025 

52 

Convergent 

4 

l-3x 

exp(-5x'") 

0 

52 

Unstable 

5 

l-3x 

exp(-5x‘) 

0 

62 

Neutrally stable 

6 

l-3x 

exp (-5x) 

0 

62 

Neutrally stable 

7 

l“3x 

exp(-5x) 

0.1 

52 ’ 


8 

l-3x 

cxp(-5x) 

0.1 

62 

Convergent 

9 

l-3x 

exp(-5x) 

1.0 

52 

poor accuracy 

10 

l-3x 

exp(-5x) 

1.0 

62 

/ 






</) = 8<X)U 
a(X) ' 3.6-3X 
e = 0.025 
DT/DX= 10 

EXACT STATIONARY 

SOLUTION (e = 0) 

O CONSISTENT B.C. 



Fif»ure 14.- Steady-state solution of a linear equation with a positive flux 
gradient and a small amount of artificial dissipation added. 
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i 


4 


i 


shows that the well -conditioned nature of the algebraic system of difference 
equations is not necessarily equivalent to the well-posedness of the differ- 
ential problem that one attempts to solve. The result of Test No. 3 is also 
indicated on Figure lA. It appears that in applying a quite erroneous 
boundary condition at x ■ 1 results a perturbation in the steady-state 
solution that is localized to a small neighborhood of this boundary. This 
is another aspect of the well-conditioned nature of this problem when 
a(x) > 0 everywhere. 

In the experiments numbered A-10, a(x) = 1 - 3x so that the sign of 
a(x) switched from positive to negative at x * 1/3. Without dissipation 
added (e * 0), ;ind a mesh of 52 grid points, the algorithm was found 
unstable even w;Lth the exact stationary solution for initial solution (Test 
No. A). However if 62 grid points are used, the i:umerlcal solution remains 
bounded (Test No. 5). This is because positive and negative elements of 
the sequence a^ * separated by a true zero in the latter case. 

However, the solution does not converge to a steady state for a different 
Initial solution (Test No. 6). If dissipation is now added (e >0), steady- 
state convergence is obtained but the accuracy of the steady-state solution 
is very poor, as shown on Figure 15 for Tests No. 8 and No. 10. One observes 
on this figure, that If e is too small, say e «= 0.1, a peak appears in 
this solution near x » 1/3. The reduction of this peak requires excessive 
amounts of artificial dissipation which degrades severely the solution 
accuracy. An energy concept can explain the existence of this peak. For 
this purpose, de.flne the following "energy" function: 

E(t) - f u(x,t)dx (8A) 





i 


. 'LL jipif ™ II' ' ^ 

ISL 




<t> = a(X)U 
a(X) « 1-3X 
DT/DX = 10 

EXACT STATIONARY 


SOLUTION (e°0) 



0 .2 .4 .6 .8 

X 


L5.- St'<?ady-state solution of a linear equation with a flux gi 
changing sign and some artificial dissipation added. 
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and consider the case e ■ 0 for which 

.1 


E' (t) ■ J u^(x,t 


)dx 


/■ 


i«|)[u„(x),x] - (Ji(u,x)}j^ dx 


= 0 (85) 

since u(x,t) is constrained to equal u„(x) at x = 0 and x -• 1. As a 
result, E(t) is a constant (nondlssipative phenomenon). Since, also, 
the characteristic curves are convergent at x = 1/3 (compression wave), 
this energy is accumulated at this point, in the limit t ^ In fact, it 
can be obtained directly from the differential equation that 

u(l/3,t) = u„„(l/3) + C exp(3t) (86) 

for some constant C, while for x ^ 1/3, u(x,t) converges to u„(x) in 
finite time. Hence this problem does not have a steady-state solution in 
the ordinary sense, unless the starting solution Uq(x) is trivially chosen 
to be U(j,(x). These experiments indicate the difficulties encountered in 
attempting to achieve the stationary solution u„,(x) by a viscosity method, 
when a(x) changes sign. 

Very similar results were obtained for nonlinear test equations. For 
these, six experiments, defined in Table 3, were conducted. Here, the start- 
ing solution u^(x) was obtained by adding to the stationary solution u„(x) 
a second-degree polynomial q(x) ■ 5x(x - 1) that is zero at the boundary 
points and negative at interior points. In this manner, negative as well 
as positive values appeared in the initial solution. Here the flux function 
<|i was chosen to depend on u only. 
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Table 3. Numerical experiments for nonlinear test equations 


Test No. 

<l>(u) 

Uq(x)-U„(x) 

£ 

Comraents 

11 

u -f u^/3 

5x(x'-l) 

0 

Neutrally stable 

12 

u + u^/3 

5x(x-1) 

0.025 

Convergent 

13 

u + u^/3 

5x(x-l) 

0.005 

Convergent and 
very accurate 

lA 

u ^/2 

5x(x-l) 

0 

Unstable 

15 

u2/2 

5x(x-l) 

0.025 

Convergent 

16 

u 2/2 

5x(x-l) 

0.005 

Convergent and 
very accurate 


For the first three cases, the wave speed a(u) = 94>/9u = 1 + > 0, 

and Instability could not be triggered. Without dissipation added (Test 
No. 11), the solution does not converge (to steady state), buf remains 
bounded. This is indicated by Figure 16 where an intermediate solution, 
obtained after lO** a.plications of the algorithm, is shown. On this figure, 
the values of u at the points Xj, Xj, Xj, . . . fall on a smooth curve, 
and so do the values of u at the points Xj, x,,, Xg, . . . but the two 
curves are distinct. This is because central space differencing does not 
couple the two subsequences (uj, U 3 , Ug, . . .) and (u 2 , u^, Ug, . . .). 

This is a known reason for requiring the use of artificial dissipation when 
a leap-frog type differencing is employed. However, if a small amount of 
dissipation is added, the numerical solution converges to a steady state. 

As on example, a very accurate solution obtained with e » 0.025 (Test 
No. 12) is shown in Figure 17. For an even smaller value of e (Test No. 33), 
exact stationary soluilon, u„(x), and mnnerical steady-state solution are 
Indistinguishable to plottable accuracy. 




0 .2 .4 .6 ,8 1 

X 


Figure 16.- An Internudiate solution of a nonlinear equation with a positive 
flux gradient and no artlflcii:! dissipation added. 
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0 .2 .4 .6 ,8 1 

X 


Flj^ure 1/,- Steady-st:ate solution of a nonlinear equation with a posit 
flux gradient and some a. Lificial dissipation added. 
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The last three experiments deal with a mcdified Burgers’ equation 
obtained by letting * u"/2. In this case, the wave speed a(u) = u 

changes sign twice at the Initial time, but not at the steady state. it no 
dissipation is added, instability occurs (Test No. 1^0* However, for r 
sufficiently large, the solution does converge. I :e 18 shows the steady- 
state solution which Is obtained for t - 0.025 (Test No. 15); again, an 
even more accurate solution can be obtained for a small value of t: , say 
e « 0.005 (Test No. 16). Here the solution accuracy Is not significantly 
degraded by the addition of artificial dissipation. This tends to indicate, 
for this pimple problen. at least , that viscosity methods converge when the 
wave speed a « does not change sig.i (at least) at the steady state. 

In conclusion, the experiments do confirm that a particular form oi 
instability can triggered when the wave-speed a(u,x) *= 3^/3u, which 
plays the role of an eigenvalue of the Jacobian matrix of a flux vector, 
changes sign at some points If central space differencing ' used at this 
point. Severe solution accuracy degradation was experienced for a case 
where the nonuniform steady-state solution was such that the alternation of 
sign in the wave- ?peed remained at the steady state. This w^s to the 
extent of making the practicability of viscosity methods questlrnablo for 
this case. However, the extension of this dramatic resiil i to the solution 
of the Euler equations is uncertain. 

C. Fmthev Comnents 

The derivation ol Section IVA above has brought a ratlon.ile, b..sed ot. 
matrlK analysis only, to explain one ot the reasons for ^he necosfltv of using 
artlflcia] dissipation. It also suggests that better stabJlltv properties 
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would perhaps result from simple modifications of the differencing at the 
points where one < f the eigenvalues of the Jacobian matrix A or B changes 
sign. 


For example, if the passage of say a^^ through zero is smooth, the 
Jacobian matrix could bo synthesized at one point from a modif 'ed eigen- 
system in which a^^ would be set equal to zero exactly. This procedure 
ensures that the matrix has purely imaginary eigenvalues only in the 

case of a scalar equation. However, a reduction in the required amount of 
adc«^'d numerical dissipation would perhaps result from applying this 
techn que. 

Another procedure, applicable to the Euler implicit method, consists 
of averaging conservative with nonconservative differencing. This gives: 


- Ax(6^ + A(S^) 

- BTrid(-(Aj_j + Aj), 0, + A^)] (87) 

which would generate a real ski"-'-r',mmetric matrix if A we.e symmetric, 
whith is rot. in general, for the Kulor equations. However, a more favor- 
able eigenvalue-spectrum can be anticipated from this. Using different 
arguments, Krelss and Oliger [19] proposed essentially this for Burgers' 
equation. Since A and E are not themselves symmetric matrices, for the 
airfoil calculation of Chapter III, using a uniform mesh in this test, the 
algorithm was found to be stable with only cwice larger time-steps when 
using this technique. 
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V. CONCLUSION 

In this work, two conjectures were first inade. First, it was claimed 
that the domain of unconditional stability of the base algorithm could be 
enlarged by introducing artificial dlssipctlou in the implicit part of the 
differencing, as well as in the explicit part. Second, it was anticipated 
that the Iterative convergence properties of the algorithm could be Improved 
by the use of larger time-steps per se, but also by the use of a cyclic 
sequence of time-steps. 

A heuristic stability analysis brought a theoretical support to the 
first conjecture. This analysis suggested that the time-step and the dissi- 
pation term added explicitly could both be arbitrary, provided the dissipa- 
tion term added Implicitly was kept sufficiently large. In particular, the 
two dissipation terms could be kept proportional to the time-step. In this 
way, the consistency condition was met, and the steady-state solution was 
independent of the time-step. This has been well confirmed by the numerical 
experiments that were conducted on a model transonic flow problem governed 
by the Euler equations. In fact, for this problem, it has nevei been pos- 
sible to find a large enough value of the time-step, for which any adjust- 
ment of the dissipation terms would not remedy stability problems. However, 
for extremely large values of the time-step, the required amount of artifi- 
cial dissipation was so large, that the iterative properties of the algo- 
rithm were degraded, although the numerical algorithm was stable. This was 
attributed to nonlinearities. For this reason, it was found that if a 
single time-step was used, this time-step could be optimized to a value 
roughly one order of magnitude larger than the one permitted by the base 
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dxi-f erencing scheme. In this manner, the modified algorithm was found to 
converge about eight times faster than the base algorithm. Even raoro rapid 
convergence was obtained by using a sequence of time-steps. With the best 
sequence, an improvement in rate of convergence by a factor of 10 (over the 
base algorithm) was observed. 

Various numerical experiments have shown that the modified algorithm 
was not very sensitive to nonoptimum parameters. In particular, approxi- 
mately the same convergence rate was obtained when using a sequence of 
eitlici. four, six, or eight time-steps. Also, the nonoptimality of the 
sequence of time-steps, for a fixed number of them, did not seem to degrade 
the convergence rate severely. 

Finally, a particular form of instability that the algorithm is sub- 
ject to has been attributed to the spaclal variation of the Jacobian 
matrices of the Euler equations. This instability found to occur when 
central space-differencing is used at a point where o. .• of the eigenvalues 
of cither one of the Jacobian matrices changes sign. Addition of a suffi- 
cient amount of artificial dissipation remedies this type of instability. 
Never thel€*.ss, two techniques have been proposed that could reduce the amount 
of required artificial dissipation. More conclusive results on this topic 
would, however, require further research. 
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VIII. APPENDIX A: MATRIX FORM OF THE FINITE-DIFFERENCE EQUATION FOR THE 

CASE OF A SCALAR DIFFERENTIAL EQUATION 

In this appendix, the (matrix) definitirn of Kronacker products and 
sums, and some of their properties are first recalled. With this back- 
ground, the finite-difference equation for the case where the implicit 
algorithm is applied to the two-dimensional first-order wave equation is 
derived in a form particularly convenient fcr the stability analysis of 
Section IIB. 


A. Some Background on Kronecker Products and Sums 
The definitions and the essential properties of Kronecker products and 
sums can be found in most books on matrix theory (e.g., [17] or [20]). The 
properties that are used in Section B of this appendix are repeated here, 
without proof, for the reader’s convenience. 

Let A and B be two square matrices of dimension JxJ and K^K 
respectively. The Kronecker product of A and B is denoted by A ^<) B and 
defined as the square matrix of dimension JKxJK given by: 


A ® B 


a^^B a^ 2 ^ 
^21® ®22® 


®J1® ^J2® 


®1J® 

^2J® 


®JJ® 


(Al) 


The Kronecker sum of A and B is defined as the matrix A ® Ij^ + Ij ® B 
where 1^^ is the mxm identity matrix. 


I 





7A 


i . 


A , 


j 


A ‘ 


The following properties are true: 

(A 0 B) ® C » A ® (B ® C) 1.A2) 

(A + A’) 0 (B + B’) = A ® B + A’ 0 B + A ® B’ + A’ 0 B' (A3) 

(A ®B)(A' ®B’) = (AA') 0 (BB') (AA) 

where A' and B' have the same dimensions as A and B, respectively. 

It follows from Equation (A4) that if A and B are nonsingular, then 
so is A 0 B and: 

(A 0 B)“l - A"^ 0 B"^ (A5) 

An important result concerning the eigeneystemt. of Kronecker products and 
sums can be stated as follows: If Aj, A 2 , . . •> Aj are the eigenvalues 
of A and Pj, ^ 2 * • • • ♦ eigenvalues of B, then the eigen- 

values of A 0 B are the numbers A^pj^ and the eigenvalues of the Kronecker 
sum of A and B are the nuniers A^ + Pj^. For both, the corresponding 
eigenvectors have the form: 

J 


^jk 


L*J 


(A6) 


where Is the mth component of the eigenvector of A associated 

to Xy and y is the eigenvector of B associated to 


t 

« f 


B. Application to the Finite-Difference Equation 
In this section, the implicit algorithm is applied to the two- 
dimensional flist-order wave equation (Equation (10)). The calculations 
art made assuming the solution u equal to zero at the boundaries, but the 
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4 


I 


result would hold for other linear boundary conditions as well. The equiva- 
lence between operator notation (Equation (3)) and matrix notation (liqua- 
tion (11)) Is expllclted. 

For this purpose, the components of the solution vector u are con- 
ventionally ordered as follows: 


U ■ ^^11* '^12* • * '^21’ '^22* • • •» '*2K’ • • *^J2’ • • •• 

^ (A7) 

where, as usual, u^j^ «= u(x^, y^.), j * 1, 2, . . ., J and k * 1, 2, . . . , K. 


J 


Let eJ and E”, E^, and E^ be the forward and backward shift operators 
for the X and y directions. More precisely: 

E^U * ^^2 1 * ^22 * • • ♦ . ^31* ^32* * * ** ^3K^ •••. 0) 

ExU - (0, 0, . . ., 0, Ujj, Uj2* . . Ujj^, . . ., * * *’ 

* ^^12* ^13* * * ** ^22* ^2 3* •••> ^J2 ^ ^J3* • • •• ^) 


E^u - (0, Ujj, 


'^21’ • • •’ “2.K-1’ 


0, u 


Jr 


'^J,K-l> 


(A8) 


vrtiere boundary conditions have been taken into account. 

For m ■ J or K, let be the mxm identity matrix, and Em be 


the following mxm matrix 


Em 


0 1 
0 1 


0 1 

OJ 


(A9) 


I 


4 


.f 

i 

\ 


6> 

V 







and Ey are 


If the matrix repreaentatlora c-f the operators E^, E“, Ey, 
denoted by the same symbols as the corresponding operators, it is i pparent 
that : 



Ej e 

(AlO) 



(All) 



(A12) 

II 


(A13) 


Clearly, these equations also hold for periodic boundary conditions, if the 
lower left comer element of Ej,, in Equation (A9) is set equal to one. 

In this case, E^ * E^^ so that forward and backward shift operators with 
the same subscript are Inverse of one another. 

Linear combinations of Equations (AlO) through (A13) and of their 
powers and use of £q[uatlons (A3) and (A4) yield the desired finite- 
difference equation (Equation (11)). 
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IX. APPENDIX B: EIGENSYSTEMS OF TRIDIAGONAL MATRICES 

In this appendix, the eigenvalues and eigenvectors of trldlagonal 
matrices are recalled. This Is done in two cases; (1) periodic boundary 
conditions, and (2) specified boundary data. 

Consider first the case where periodicity is enforced at the boundaries, 

so that the general tridiagonal JxJ matrix has the following form: 

’be 
a b c 
a b c 


A = Trld(a,b,c) * 


Lc 


a b c 
a b. 


(Bl) 


Define a sequence of vectors Xj (j = 1, 2, . . . , J) by 

X^j - l//j exp(mi6j) (B2) 

where 6j = 2n(j - 1)/J and m = 1, 2, . , J. If Equation (B2) is also 
applied for m » 0 and m ■ J + 1. the periodicity boundary conditions: 


X , = X^ , 

0, j J,j 

Y *- X 

1. j J+i,j 


(B3) 


are found automatically satisfied by X j . Now compute AX j : 


(AX, ) - A , X, , 

J m mk kj 




(B4) 


where 


Aj " a exp(-10j) + b + c exp (16^) 


(B5) 


I 





I 


this shows that is an eigenvector of A associated to the eigenvalue 

Now compute the following Inner product: 




where a 


2(r(k - j)/J. This gives 


<Xj,Xj> - 1 


for all j, and for j / k: 


> - i expdc.) 

y \ J ^ expCicx; - 1 


Hence, the eigenvectors (j * 1, 2, . . . , J) are orthonormal* The matrix 
X that diagonalize^ which contains tnese eigenvectors for column vec- 


tors, is thus unitary: 


X* (adjoint of X) ^ X ' 


Consider now the case where the components of the solution vector are 
const ralm'd to be zero at the boundaries* For this case, the general J<J 
tridiagonal matrix is given by Equation (Bl) In which the upper-right and 
lowc’-r-left corner eleraents of the matrix on the right-hand side of this 


equation are set equal to zero* Then, define a sequence of vecLors X. 


(j - 1, 2, * . *, J) by; 


1) (/a/c)™ sin m6 , 


(BIO) 


where 0^ =» + 1) and m =» 1, 2, • . . , J* If Equation (BIO) Is also 

applied for m - 0 and m ■ J + 1, the boundary conditions 


^ 4 * 0 

Ofj 

are found automatical] y satisfied by * Now compute AX^ : 


(BID 





\ 




aX . + bX . + cX . , . 
m-l,j inj m+1,3 


j mj 


sin(m - l)0j sin(m + l)0j 

"j " ^ sin m0j (3^3. 

- b + y^v. cos 0j 

which is found independent of m as expected. This shows that Xj fs an 
eigenve -or of A associated to the eigenvalue X j . 

Nov; consider the particular case where a, b, and c are real, with 
also a » c, and redefine the eigenvectors of A to oe Clearly, 

those are (real) orthogonal since A is real symmetric in this case, so 


that: 


En plus: 


= 0 (for j ^ k) 


(«rS) "mZl 


(1 - cos 2mep 


J - R 
J + 1 


where 


R = Re(S) 


S » exp(m-*a) 
m* 1 
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in which a 


2M/(J + 1). 

S - 


Computing S gives: 


cxp(ia) 

exp(ia) 


J-l 

cXp(miu) 

m*o 

exp (Jig) - 1 
exp(ia) - 1 


4 Jet 
sin Y 

*tpi(J + l)^ — 

2 . a 

sin *2 


(B17) 


so that: 

cos (J + 1 ) ^ sin ^ 

R ^ (B18) 

sin *2 

Also : 

C03(J + 1) ^ - cos ~j - ( ')^ 

sin 4^ • ^In ^ sin 4 

so that R “1 and *^1. As a result, the matrix ^ that diagon- 

alizes A and contains the eigenvectors ^ for column-vectors is 
erthogona] : 

= I 


This matrix Is pu plur syranetric : 


■; - f (B20) 

In particular, ' diagonalizes the smoothing operator Dj^ of Equation 0 2) 

“ *:(2I + 2K)£; - 2(1 + CKO (B21) 


where iC * Dlag(cj) aid Cj ■ cos 

Now consider the matrix * Trld(-1, 0, 1), in Eq>iatlon (12). In 
view of Equation (T.IO), it appears that this matrix l-s diagonalized by 
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^ t 


X - nr. (P22) 

where D is th. JlaKoual mitrix with mth flgenvalue fqiwl to 1~. Hotui . 

X-» - r’n--‘ - -.D - -*D* » X* (B 23 ) 

showing that X Is thvn unitarv as expected, since Is real skew- 

syinDetrlc. The eigenvalues of C . are given hv Fquatlon (Bni. and this 
permits us to write t ^ In the f.>l lowing form: 


' X(1K)X 


- 1 


wheie the diagonal matrix K Is d«*flned in Kquatlo:: 
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X. APPENDIX C: STABILITY CONDITION FOR PERIODIC BOUNDARY CONDITIONS 

and second-order smoothing 

In this appendix, Equation (39) which expresses the stability condition 
for the case of periodic boundary conditions and second-order smoothing, is 
derived. For this purpose, it is recalled that the satisfaction of Equa- 
tions (33) and (35) constitutes the necessary and sufficient condition for 
stability. 

For the case considered. 


dj - dj ■ 2(1 - cos Oj) 
dj^ ■ d^ = 2(1 - cos ^j^) 


(Cl) 


where 


2'(j - 1)/J (j - 1. 2 1) and 


(k « 1, 2 K), and it is convenient to let: 


- :-(k - n/K 




dj - 4f 


‘*k * ‘‘k ■ 


/A 

“/A 


1 


(C2) 


In this tnunner, the new variables \ and n, which are not subscripted for 
notatlonal slmp)iclt7, take their values in the interval (0,J]. Kqua- 
tion (3h) then becomes: 

g , ( 0 ) * ^ (f + n)‘^ - 2‘’f(r, + n)f0r(f. + n) + 2] 

j ** 

+ 2c(e + n) + Om - m)*' 

• o" [(r - + n)^ - + n)l (C3) 


f ! i ! I 1 I 
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where u ■= (6 - l/2)/('^. As a result. Equation (35) becomes: 

r(f:,n) < ue (C4) 

where 

/- _ -\2 . 

(C5) 


t(t.n) - (f. + n) - ^ 


(C7) 


The satisfaction of Equation (C4) for all feasible values of f, and n is 
equivalent to the condition: 

S(r,f) < (C6) 

where 

S(c,f) = Sup f(f.,ri) 

0 i f, , > i 1 

i 

Somr sltnple conclusions can be drawn at first. For this, note that. If 
c • r , f ,n) < 0 and Equation (C6) is satisfied. Thus, the algorithm is 
unconditionally stabli- for every value of 0 (provided Equation (33) holds) 
for the case where the equat“*on 

IV - + bu^ « 

is differenced In a time-accurate manner. Also observe that 
f(l,0) ■ (c - c)‘/4 2 0* tlvat, unless perhaps this value is zero, S(t ,r) 
is strictly positive. Consequently, for trapezoidal timo-dlfferencing 
(6 •= 1/2), since \i • 0, letting * ^ 7 constitutes the only way of 
enforcing unconditional stability ^see Equation (C6)). For tbls reason, 

0 > 1/2 and m > 0 are assumed in vhe remaining. 

The next step consists of the determination of S(t,r), that is, the 
maximization of f(C,i) for fixed values of i and r. Inspection of 
Equation (C5) Indicates that f(^,n) is a homogeneous function of ^ and n 


E 
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\ 


i 


of degree one. Applications of Euler’s theorem for homogeneous functions 
gives the following identity: 

f«,n) - II C + |i n (C9) 

This shows that If there exists a local maximum of f(C,n) at a point of the 
open square ]0, 1[ x ]o, 1[, this maximum is equal to zero. In view of 
Equation (C4), it appears that such maximum, if it exists, does not yield 
any binding condition for stability. Hence, the maximization of f(C,n) 
can be reduced to the one over the boundaries of the square. Since, en plus, 
f(C»n) is symmetric in C and ti, only two of these boundaries need to be 
considered. These are the segments: (1) 0 < C < 1 and n - 0, and 

(2) C “ 1 and 0 < n < 1. 

For n ■ 0, 

f(c,0) - (CIO) 

It is maximum at the point C “ 1, which also belongs to the second se^ 'sent. 
Consequent ly, 

S(c,c) » Sup 4>(ri) I 

} (Cll) 

0 < n < 1 1 

where 


(Kn) - f(l,n) 

• (^T^) + 1) + 77^ - (C12) 

Differentiating Equation (C12) gives: 

* (^n + 1)) ” “"1 <C13) 

where oj »= 2e/|? If o)il or uj<2. <})’(n) does not change sign 

when n increases from 0 to 1, and (j>(n) is maximum at either one of the 


/ 


,s\ ' i y 


A 



. i 
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two endpoints If 1 < u) < 2, (n) < 0 for 0 < n < w 1 and > 0 

for w - 1 < n i 1, so that <j)(n) is again maximum at either one of :he two 
endpoints . Finally : 

S(e.e) - Max[({(0), 0(1)] (C14) 

One obtains: 


> 0 


♦ <o, . (L^) 

((>(1) “ c - e) 

Equation (C6) then breaks into the following three Inequalities; 


"2 ■ V* i E 


(C15) 


(C16) 


Using the definitions of r and c given in Equation (C2) and making a few 
simpllflcntlons yields Equation (39) (Q.E.D.). 

Remar k : In this derivation, .10 case has been made of the "hape of the 


function h(0) for which d 


1 




Equation (C16) applies to all the cases where the same type of smoothing is 
applied implicitly and explicitly, provided r and e are defined by: 


cr . / X 

1 max 


E = r /X 

® max 


where it is assumed that the eigenvalues cf the smoothing operator vary from 
0 to In particjular, for fourth-order smoothing applied explicitly 

and implicitly one wotild let: 


16 H 


A 

16 

1 

16 


and apply Equation (C16). 



I 




86 


XI. APPENDIX D: STABILITY CONDITION FOR PERIODIC BOUNDARY CONDITIONS 

AND FOURTH-ORDER SMOOTHING 


In this appendix, Equation (40), tdiich expresses the stability condi- 
tion for the case of periodic boundary conditions and fourth-order smoothing, 
is derived. For this purpose, it is r<!called that the satisfaction of 
EquJitions (33) and (35) constitutes th«! necessary and sufficient condition 
for stability. 

For the case considered. 


dj - 2(1 - cos 0j) , 
2(1 - cos 0. ) , 


dj ^ d^^ 


‘^k 


f 


(Dl) 


't k^ 

r 

where 6j = 2it(j - 1/J (J « 1, 2, . . ., J) and •» 27i(k - 1)/K 
(k •= 1, 2, . . . , K) , and it is convenient to let: 

^k 


4f . 

d^ = 16£;2 

4n , 

d; = 16n' 
k 

0F / ' 


r/16 



(D2) 


In (his immner, the nc^w variables r, and n, which are not subscripted for 
notationa] simplicity, take their values in the interval [0, 1], Equa- 
tion (36) then becomes; 


g^k^®^ - o2?-'(c2 + n')^ ~ 2er(5? + t|2)(0e(C + n) + 2] 
+ ?F(C'^ + n‘ ) + 0 ^r^(^ - p)-’ 


Def 1 ne: 


g(C.n) 


2 


(D3) 


(D4) 
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and 


* g2(5,„) _ ^^(5? + v2) - (D!,) 

where u *= (0 - 1/2) /f'^, so that 


gjj^(e) - 4e2f(e,n) (d6) 

As a resuJt, Equation (35) becomes: 

f(C,n) < 0 (D7) 

The satisfaction of Equation (D7) for all feasible values of C and n is 
equivalent to the condition: 


S(e,t) < 0 


(D8) 


where : 


S(e,c) =» Sup f(^,n) 

0 < 5, n < 1 


(DP) 


Thus, the problem consists of maximizing f(?,n) over the closed square 
[0, 1] X [0, 1]. For this purpose, one first looks for stationary points 
of f(C*ri) that belong to the open square. At these points, if any exists. 


— « 2g(?,n) ||- - 2gtC - . q 


Tn ’ 


•r^ — 3urn — f2 


(DIO) 


- 0 ; 


so that, in particular 


0 - ^ 

.>n 3C 


- ;'g(f.,n)0^ ■ If ) ■^ - e 2)(C - .1) (Dll) 

- (f, - n)[-2eg(f;,n) + 2vie - e^] 

which requires the satisfaction of at least one of the following two 
eqwitlons : 



38 


5 * n 

2vit - 2eg(^,n) » 

If Equation (DIO) holds, It is also true that: 


0 


iL r M 

3n ~ ' 3£; 


- 2g(C,n)^n T^' - C + 2ye(r,‘ - n^) 

- 2g(e,r))[c(n^ - - I (r, - O] + 2u£(f;^ - n^) 


(D12) 

(D13) 


■ ii' - n^)(2ye - 2cg({;,n;] + cg(5,n)(£; - n) (D14) 

Now, assume that f(^,n) Is stationary at a point (C,n) of the open square 
that does not belong to the line « « n. Then, at tils point. Equa- 
tions (D13) and (D14) must hold simultaneously, so that: 


0 e(C + n) + g(f.,n) 


\ [r(^- 


+ n^) + c(f. + n)l 


(015) 


If the trivial case e •« F * 0 is eliminated, the satisfaction of 
Equation (D15) ls Impossible tor f, 0 and n > 0. This brings a contra- 
diction to the assumption r ^ n. Consequently, if f(f,T)) is stationary 
at a point (r. of the open square, this point belongs to the Line C * n. 
Hence, it suffi*'es to enforce Equation (D7) on the boundaries of the square 
and Its diagonal segment 0 < i » i) < 1. Observing that is sym- 

metric in r, ami m permits us to further reduce its maximization to the 
one over the following three segments; (1) 0 < r, * n < I, (2) 0 f. ^ 1, 
n ■ 0, and (3) . ■ 1, 0 < n 1, For this reason, three cares will now be 


examined separately . 
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Case 1: 0 < 4 ■ n < 1 


Define 


^(4) - f(4.4)/42 

■ (e4 ~ - (2\ic + c^) 


(D16) 


SO that 

= 2e(cK - e) (1)17) 

When i increases from 0 to r/e, 4>(0 decreases, and since 4>(0) ■ -2 vje < 0 
remains negative. For values of f, greater than e/e, is an increas- 

ing function of f;. The satisfaction of Equation (D7) over the considered 
segment is hence equivalent to the condition 

<t(l) < 0 (D18) 

This gives the following necessary condition for stability: 

C > ~ (c - 2v) (D19) 

Case 2 : 0 i 4 < 1, n * 0 

For this case. Equation (i'>7) takes the following particular form* 

- UC4-' < C (D20) 

which is equivalent to: 

(cC - r)‘^ < (D21) 


or 


?v^, - < c < €i + 2»^uf (D22) 

The binding case for the inequalltv on the left corresponds to C » 1, For 
the inequality on t!ie right, the binding case corresponds to r, * 0. This 
gives the following necessary condition for stability: 

€ -• 2/ImF i e < 2^^ (l)23) 




tiM 'aw 
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Case 3; * 1, 0 < n < 1 


For this case, one defines 


<J'(n) * f(l n) 


i I?(n^ + 1) - t:(n + 1)]^ - ;.c(>r + 1) - e^n 


(D2A) 


which Is a fourth-degree polynomial In n. The satisfaction of Equation (D7) 
over the considered segment Is equivalent to the condition 

<t'(n) < 0 (D25) 
for 0 < n < 1. It Is assumed that Equations (D19) and (D23) hold, so that 
iJ^(O) and v(l) are both nonpositive. Hence, if an additional condition must 
be enforced, it must be of the form: 


4’(n ) < 0 


(D26) 


where 0 - n* < 1, and 


- 0 


(D27) 


(foi a local maximum). Dividing the polynomial 4 (n) by its derivative 
tj^'(ii), according to d»‘creaslng powers of ii, produces a quadratic remainder, 
q(n). This gives the following identity: 


li'(n) * (an + b)v'(t!) + q(i;) 


(D28) 


where 


q(n) = an’ + ibi + > 


(D29) 


The calculation of the coefficients a, b, a, d, and > gives: 

a = 1/4 


/(8F) 


(D301 


a ^ |4f(r - - 2p) - f’| 


■ 8cr| 

> = If (c - t ) ■ - 4uc 1 - - ■ 


(DID 
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li 


i 



where terns between blackens [ ] are nonpositive, as a consequence of 
Equation (D23). Hence, It is directly apparent that 3 < 0 and Y 0, 

This implies that q(0) - Y < 0. Now compute q(l): 
q(l)=a+3 + Y 

(D32) 

< - 4 et; - Spe) + (-4?^ - See) + 7 (e - e)^ - pc 

16 16e A 

where some nonpositive terms have been neglected. After a few simplifica- 
tions, one obtains: 

q(l) < f (i - 2e) - I ue (D33) 

Thih indicates that in case c < 2c, q(l) < 0. If now, c - 2e > 0, the 
satisfaction of Equation (D19) requires that e - 2e < 2u, so that 
q(l) i “Ur/2 < 0. Thus, in all cases, q(l) < 0. Now, if a < 0, since 
3 < 0 and Y £ 0, q(Ti) < C for 0 < n. If a > 0, q(ri) achieves a mini- 
mum at the point 

n = - > 0 (D34) 

Whet h r u belongs to the interval [0, 1] or not, since q(0) and q(l) are 
both nonpositive, so is ,(c) for all values of t) in this interval. 

Hence, in both cases, q(n) <0 fi. : 0 < n 1. 

Hence, if it happens that il'(n) admits a local maximum at a point n* 
of the open interval 10, 1[, tlien, as a consequence c” Equations (D27) 
and (D28), the following is true: 

= q(n*) < 0 (D35) 

This shows that < 0 for all values of ii in the open Irterval ]0, 1[, 

pro\'lded l.quatlons (DJ9) and (iJ23) hold (sufficiency). 
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In conclLSion, the satisfaction of Equations (D19) and (D23) consti- 
tutes the necessary and sulfic^ent zondltlon for thf* unconditional stability 
of the algorithm for the case considered. Replacing, In these equations, 
z and G by their definitions, given in Equation (D2), yields the desired 
equation (Equation (40)), 

Reim irk 1; 

In this derivation, no case has been made of the shape of the function 

h(0) for which d. « h(6.), d! » h*-(0 ), d. - h(0, ), d* = h^(0 ). Hence, 

J J.1 kK k 

Equations (D19) and (D23) apply to all the cases where the explicit smooth- 
ing operator (or Uy) is the square of the implicit smoothing operator 
(or Dy)j provided c and c defined by 

(D36) 

where it Is assumed that the eigenvalues of (or D^) vary from 0 to 
^max* However, this has apparently no application, since the next step 
aftt'r tne combination second-order implicit/fourth-order explicit smoothing 
would be the combination fourth-order implicit /sixteenth-order explicit 
smoothing, which is intpractical . 

Remirk 2 ; 

Note tnat Equations (D19) and (D23) are equivalent to Equation (C16). 
This is because the same boundary values of f, and n give rise to the 
stability conditions fo^* fourth-order smoothing as for second-order srKkoth- 
ing. Hence the two cases only diff^'r by scale factors which have been elim- 
inated from Equations (C16), (D19), and (D23) by a ^ropriate definitions 
the parrmeters e and e. 


e = 0G /X \ 

1 max I 

e “ I 

® max / 



} 
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XII. APPENDIX E: EFFECTIVE EIGENVALUES OF THE SNOOTHING OPPlATOixS 

AT LARGE CCURAI^’T Iv/MBERS 


The *'etfectlve eigenvalues/* and d j , of the smoothing operators, 
and or were Introduced in the development ( l the stability 

analysis for the ca>,^ of sped- led boundary data and large Courant numbers 
(Section IIB3). These effective eigenvalues ait evaluated in this appendix. 
For this purpoi i, on recalls that: 

(El) 

dj - (E 

where the various matiices are of dimensions JxJ and defined bv; 


Trid(-l,2,-l) 


D; - Dx or Dx' 


X^j « /2/(J + 1) i® sin mO^ 


I 


x;j = v/2/(J + 1) (-l)J sin je 

where * jir/(J + 1) a.id m, j = L, 2, . . . , J. 
Appl'^lng Equation (Bj 3) to the case where: 




2 

1 - c 


(E3) 


(E4) 


gives 


2r + i'-(l - t'‘’)c?.a 6, 


■ i cos + 2f + Olf‘) 


(E5) 


This shows that when the matrl >; - Trid(-l,0,r Is perturbed by raatvlx 

IS a small parameter, the first-oiuer 


ePy ■ rrld(-c, 2f, -e), where 



i f f J I i -1 ' \ - \ t ' < ■ 1 ' ♦ t • 


\ i !■ • •1 i > ■- -4--f ..,-4 
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perturbations brought to the eigenvalues of have the coefficient 2. 

This is equivalent to saying that: 


dj .2 


(E6) 


as one could compute by expliciting Equations (El). But according to Equa- 
tions (B21) and (B22); 

S, - X-‘D^ 

' CD (21 + 2CK5)DC 

= 21 + 2U*KU (E7) 

where C, D, and K are defined in Appendix B, and 

(E8) 

is another unitary matrix. One concludes that the diagonal elements of 
U*KU are all equal to zero. Considering now the case where D' = D ^ and 
squaring Equation (E7) yields: 


K - X-1d^2x 


Hence, for this case; 


where : 


41 + 8U*KU + 4U*r2u 


dj . 4(1 




(E9) 


(ElO) 


(Ell) 


m' 


in which c^ - cos 6j„. It has not been possible tc evaluate wjj explic- 


itly. However, it appears from Equation (El.l) that wjj. and consequently 


f . r T ■ 


1 ? - ' t ^ - 1 < • i 1 = '■**'’4 
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t 

^1 


dj is teal and positive. Moreover, since U Is unitary and symmetric, the 
following is true for every value of j: 

[2 - 1 (E12) 

Since IcjJ < 1, the following bound holds: 

“ ^ “jj ^ 
and consequently: 


A < dj < 8 


(E14) 


The remainder of this appendix shows that the maximum value of dj does 


converge to 8 in the limit of a mesh refinement. For this, let 


w = Max w 

J 


jj 


to make the claim equivalent to the follcwlng statement: 

lim w ■ 1 


(E15) 


(E16) 


The simplifying assumption that J is odd is made, and one lets 

J + 1 = 2v (E17) 

so that * rrv/(J + 1) ■ tt/ 2. Since w < 1 (see Equation (E13)), it is 
sufficient to show that: 

11m Wyy ■ 1 (E1&) 

J->« 

To evaluate one first computes using Equation (E8)r 

V - 

J 


2 1 

7~T~i ' ? J ^ *"®J si** J 

j-l 


e. 


(E19) 


* 4 


» 4 


* 4 t ( 1 


.-I 
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where the values of f and D were taken from Appendix B. Recall that 
0y * tt/ 2 to simplify Equation (E19) as follows: 


U ■ — ^ — S j^2k-l 
^mv j + 1 ^ 


sln(m02j^_j) • (-1) 


k-1 


k-1 


i sin[(2k - 1)0^] 
k-1 


In, 

J + 1 ^ 


where 


Im “ 


c^m “ ]C exp [1 (2k - 1)6^] 
k-1 


One first computes Og, as follows: 


v-l 


» exp(10g,) exp[q(210g,)] 

q-0 

exp(2vi6g,) - 1 
* exp (2iegj) -1 


exp(2vl0g,) - 1 


21 sin 0 


m 


Note that 2v0gj - ffm, so that 


. 1 


ra 21 sin 0 


m 


r^L±J=H 

^ 2 sin 0 


m+1 


m 


(E20) 


(E21) 


(E22) 


(E23) 


(E24) 
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and 


i ■ 1 + -C-U 


afl 


^mv j + 1 sin e 


m 


(E25) 


One now applies Equation (Ell) to the case j * v to get: 


"vv “ E l^mvl^ ® 

ofI 


m 


E KX - L lUmvP 6m 

ta»l m»l 


- 1 — ^ 


- 1 - 


j + 1 


(E25) 


where Equations (E12) and (E17) have been used. Equation (E26) validates 
the statement made In Equation (E18) (Q.E.D.). 
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XIII. APPENDIX F; STABILITY CONDITION FOR SPECIFIED BOUNDARY DATA 

>«) SMALL COUIUWT NUMBERS 


The stability analysis developed in Section IIB3, has shown that when 
the data are specified at the boundaries, the following approximate stabil- 
ity condition applies to the case where the Courant numbers and are 


small: 


where 


2(1 + e^d^Kl + e^d^) - e^Cdj + d^ > 0 


dj * 2(1 + cos 9j) 

d^ ■ 2(1 + cos 9j^) 

in which 6^ ■ irj/(J+l) (J - 1, 2, . . . , J) and 0j^ ■ Trk/(X + 1) 
(k « 1, 2 K), while 

“j ■ 


<n) 


(F2) 




■“k " 


(F3) 


depending on whether second-order or fourth-order smoothing is applied. In 
this appendix, the condition expressed in Equation (FI) is expllclted for 
these two cases. For this, one lets 




4C 

4n 

e/4 


1 


(F4) 


and 


dj - 4? or 16^2 
d^ ■ 4n or 16n^ 
e- - F/4 or e/16 


(F5) 
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In this manner, the new variables i and n, which are not subscripted for 
notatlonal slopllclty, take their values In the interval fO,l], and Lqua- 
tlon (FI) beccnes: 

2(1 + eO(l + en) - > 0 (F6) 
where p Is the order of the smoothing applied (p ■ 2 or 4). 


A. Second>*Order Smoothing 

In the particular case where second**order smoothing is applied 
Equation (F6) becomes, after rearrangement: 

I -2e <2 

C + n 

Hence, we are led to determine the function: 

f(e) - 


Inf g(C,n) 1 

0 < e, n s 1 I 


where 


g«.n) - 2 


- C + n I 

- n - C I 


(P ■ 2), 


(FV) 


(F8) 


and to write the stability conditions as follows: 

? - 2e < f(e) 

For this, one lets 

X 

y 

In this way, the domain of study (see Figure 19) is now defined by: 

0 < X < 2 

|y| £ y„,<x) 


(F9) 


(FlO) 


(Fll) 


(F12) 


where 2 - x depending on whether x i 1 or x i 1. 

Computing 
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*(s.n) . h(x,y) - (ns) 

Indicates that for given x, h(x,y) is minimum at y ■ 

4>(x) « h[x,y (x)]. For 0 < x < 1, y^.„(x) ■ x so that (Ji(x) - 2/x 
toax ukax 

which achieves a minimum value of 2 at x * 1. For x > 1, y„„(x) * 2 - x • 

max 

so that: 


♦ (x) 


4 -f c^[x^ - (2 - x)^l 
2x 


(F14) 


Hence, if e > 1, <p(x) is nondecreasing over the Interval [1,2], so that 
Hx) is actually minimum at x ■ 1. This gives 

f(e) - 2 for e > 1 (F15) 
If now. Instead, c < 1, <|i(x) decreases when x Increases from 1 to 2, and 
i^(x) achieves its minimum at x > 2, so that: 


for stability: 


f(e) - 1 + c2 

for 

e < 1 

(F16) 

results with Equation 

(FIO) 

yields the following 

conditions 

e > i/T ^ 1 

C - 1 

for 

for 

e < 4 
e i 4 

(F17) 


Replacing e and e by their expressions in terms of and given in 
Equations (FA) and (F5), yields the desired equation (Equation (56)). 
Finally, the remark that was made on Equation (C16) in Chapter IX also 
applied to Equation (F16). 



\ 

! 


B. Fourth-Order Smoothing 

In the particular case where fourth-order smoothing Is applied (p A), 
Equation (F6) becomes, after rearrangement: 


!(1 + tl 




(F18) 


Hence, If one defines, for this case, a function g(C,n) by 


g(C.n) 


+ n2 


(F19) 


and a function f(e) by Equation (F8) the stability condition takes a form 
similar to Equation (FIO), which Is 


e < f(e) 


(F20) 


New variables x and y are also defined as In Equation (Fll), and the 


domain of study Is still given by Equation (F12). Here, 


8(5, n) - h(x,y) 


4 + e^(x^ - v^) -f 

+ y^ 


(F21) 


Note that for given x, ly| < x In the domain, so that the numerator of 
h(x,y) Is positive and decreases when |y| Increases. The denominator Is 
also positive but It Increases with |y|. Thus again: 


f(e) = Inf (|i(x) 

0 < X < 2 


where 


(F22) 


«Kx) - 


(F23) 


For 0 < X < 1, yj^jj(x) • X, dO that: 




(F24) 


which achieves a minimum value of 2(1 + e) at x ■ 1. For 1 < x < 2, 


yjjjjjjj(x) ■ 2 - X, so that: 
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_ 4 4- £^[x^ - (2 - x)^] + Acx 
x2 + (2 - x)^ 

. 2 +, e),x t ,l_-^5. (F25) 

x2 - 2x + 2 


and 

(» ■ ! ^»(x) . (e2 + e)(x2 - 2x + 2) - 2(x - DfCe^ + e)x + 1 - e2] 

- -(e2 + c)x2 + 2 (c 2 - l)x + 2(1 + e) (F26) 


Note that <ti'(x) has two real zeros given by: 

- 1 ± >^(e2 - 1)2 + 2e(l + e)2 


Xj.Xj 


e2 + e 


(F27) 


Clearly Xj < 0 and X2 > 0, and (|i(x) changes sign at most one time (at 
X2) In the Interval [1,2]. Since, also, <Kl) ■ 2£(e + 1) > 0, <('(x) 
achieves Its minimum at either x 1 or x - 2. Computing 

^il) - 2(e + 1) 

4.(2) - (e + 1)2 

Indicates that 41 (1) < 4.(2) when c > 1. In view of Equations (F20) 
and (F22), the stability condition Is written as follows: 


(F28) 


e > - 1 

for 

e < 4 1 

e > 1 - 1 

for 

e > 4 i 


(F29) 


Replacing c and e by their expressions In terms of and given In 
Equations (F4) and (F5), yields the desired equation (Equation (57)). 

The remark that was made on Equations (D19) and (D23) In Chapter X 
(Remark 1) also applied to Equation (F29). Also note that Equation (F29) 
Is identical to Equation (F16) for a reason given In Chapter X (Remark 2). 



I 
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XIV. APPENDIX G: ON THE EIGEIWALUES OF THE MATRIX T 

The stability analysis of Section IIB was developed for th.-. aje i a 
scalar, linear partial-differential equation with constant coefficients. In 
this way, the convective derivative operator 0,^, then (real) 8kew~f»-r;anecric , 
had purely Imaginary eigenvalues, and this property was crucial to the 
derivation. In Chapter IV It was shown that for the Euler equations, this 
property Is most likely to be true when the following matrix 


‘0 a 
«1 


2 

0 a- 


S2 0 


*J-2 


a 


‘J-l 


<G1) 


has real eigenvalues. (In Equation (Gl), a^ represents the value at 
“ (J ” l)Ax (J " 1, 2, . . ., J) of any one of the four eigenvalues of 
the Jacobian matrix A (for 2-D flows), and Fj Is subscripted to Indicate 
Its dimension.) In this appendix, the question of whether T. has real 
eigenvalues Is Investigated. 

Three cases will be examined successively. They are; Case I — 

Sj, Sj, . . •* Oj all nonsero and of the same sign; Case II — there Is 

at least one true zero In the sequence at every alternation of sign; and 

Case III - there Is at least one value of J for which 0* 

All three cases could be studied by analyzing the characteristic poly- 
nomial of r^, that Is 


Aj(X) - det(Fj - XIj) 


(G2) 





I«IMNI| MM' ' 
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where I. is the JxJ identity matrix, but for Cases I and II simpler 
arguments, which yield the desired Information, are preferred here. 


For Case I, define u sequence 5j by; 


V> ■ ')Vv^ 

Clearly the sequence 6^ Is well defined, 
so that the diagonal matrix 


(G3) 

2, . . J - 1) 

and none of its elements is zero, 


7 - Diag(6j) (GA) 

Is nonsingular, and 

7-1 «= Diag(6JO (G5) 


Then perform on the matrix Fj 


7“ir,7 


the following similarity transformation: 


- Dlag(6~l)Trid(aj_j,0,aj^^)Diag(6j) 


(G6) 

(G6) 


The matrix F Is found real symmetric. As so, It Is diagonalized by an 

J 

orthogonal transformation 0 and has real eigenvalues. These eigenvalues 
are also those of the matrix Fj which Is similar to the matrix F^. This 
case Is hence favorable to the stability of the Implicit algorithm. (More 
Information about the eigenvalues of the matrix F^, In Case I, Is gi^'en In 
the remarks at the end of this appendix.) 

Using a continuity argument. Case II can be treated as an extension of 
Case I. For this, an undetermined scalar parameter x Is substituted in 
place of every zero appearing in the original sequence a^ , all the other 


t "**1 1 »•» ^(f:^ "BaCTl “iW(< swiMf •=SW <)• j s 
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elements of this sesuencc L«.ing maiptained the samp. Tn this way, the 
matrix F becomes a I'lrst-dfgree polynomial of x whose elgenvaltes for 

J 

X • 0 are required. Th,’. coefficients of the charac*-erlstlc pclynomlal 
Aj(X) are simply some polynomials of x, that is, continuous functions of 
X, so that: 


Aj(A) 


- 11m A (X) 

x*o x^o 


(C7) 


Put for X j* 0, Equations (G3) through (G6) apply. Hence, for x ■ 0, the 
matrix fj (here defined by the last line In Equation (G6)) and the matrix 
Fj may not be .similar, but have the same characteristic polynomial, and 
thus the same eigenvalues. To analyze these eigenvalues, suppose, for 
example, that ■> 0 for some r. In the neighborhood of this element. 


the matrix f.. has the following structure: 


O -/a ~a T 

T r-2 r-l 

O 


r-.?®r-i 


O 


o o 


o 


o 


Vw 


r+2 




rfl^r+2 


O 


Here, the Mtrlx 1b found block-diagonal. Any of these blocks is a 

real symmetric (trldlagonal) matrix (because consecutive nonzero elements 
of the sequence a^ a^e Df the same sign), and separated from the next one 
by at least one row and one column of zeros. In view of this, let be 

K 



1 } ■ 1 ' 




the ^orthogonal) matrix that diagonalizes the kth block, and fi be the 
(orthogonal) block-diagonal matrix where the blocks are taken to be 
Imj» ^®2’ ^3’ ^ Identity matrix and mj^ 

the number of zeros that separate posicive from negative elements of the 
sequence a^ at the kth alternation or sign. Clearly the matrix 0 
diagonalizes the matrix f whose eigenvalues r.re then found to be those 
of the blocks in the "d agonal" of the matrix taken together, with 

en plus the eigenvalue 0 added Em]^ times. These eig**nvalues, which are 
also those of the matrix Pj, are all real as is desirable for stability. 

For Case III the characteristic polynomial ^j(^) of the matrix Pj 
nfc..i8 to be analyzed. For this analysis, expand the determinant of the 


matrix 




^J+1 ■ ^^J+1 


along its last column to get: 


^j-i - '"j-i 




and finally: 


“ -AAj(X) - aj3j^jAj_j (X) 


From this result, one can derive the following formulas: 


J 


J 

.4 

► ‘.i 
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where n is a natural lnt«*g«r> and Pg^Cx) and Qg,(x) are polynomials ot 
degree m, with leading term x™, and which satisfy the following recurrence 
formulas: 

To see this, use induction. Compute 

l-X a, 


Ajd) - 


■1 

-X 


- ajS2 


and 


A,(X) 


“X 0 

a . *“X a . 


0 Sj -X 


« X(x2 - Sjaj) - ajC-ajX) 
■ -XIX^ - (a^aj + aja,)] 


Clearly, if one defines: 


Ao(X) - 1 


J'o' 


(G12) 


(G13) 


(G14) 


,(x) • Qo(x) • 1 
Pj(x) - X - aja 2 
Qj(x) - X - (»i*2 * *2®3^ 

Eqwitlon (GIO) holds for m - 0 and m • 1, .snd Equation (Gll) holds for 
m - 0. Now suppose that Equation (GIO) hclds for m - r for some r > 1, 
and Equation (Gll) holds for m • r - 1, Then applying Equation (G9) with 
J • 2r + 1 gives: 





•( V f • "I - "I i>^.l ' f 'i ' ' i ^ 
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“ ®2rfi®2r+?^2r^^^ 

- x2q^(x2) - (G15) 

Hence, Equation (GlOa) holds for m - r + 1 provided P^^(x) Is defined 
according to Equation (Gila) In vhlch n Is set equal to r. It also 
appears from this equation, that since P,.(x) and Q^(x) are polynomials of 
degree r, with leading term x’^, (*) a polynomial of degree r + 1, 

with leading term x . Similarly, 


^2r+3^^^ " ”^^2rf2^^^ “ ®2r+2®2r+3^2t+/^^ 

--X[P^j(x2) + a^^2a^^3Q^(x2)] 


(G16) 


Hence, Equation (GlOb) holds for m ■ r + 1 provided P^j(x) Is defined 
according to Equation (Glib) In uhlch m is set equal to r. It also 
appears from this equation that since P . , (x) and Q (x) are polynomials of 

ITtI X 

r 

degrees r + 1 and r, respectively, with leading terms x and x , 

respectively, la a polynomial of degree r -t 1 with leading term 

^ l- 2 

X . This shows that Equation (GIO) holds for m <■ r + 1 and Equation (Gll) 
holds for m ■ r; hence, they both hold for evfery value of m. (Q.E.D.) 

In view of Equation (GIO), the following two conclusions can be drawn: 

(1) If X Is an eigenvalue of the matrix Tj, then -X Is also an eigen- 
value; and (2) the eigenvalues of the matrix ij are all real If and only 
If the roots of Pm(x)» m case J ■ 2m, or Qj,(*)» In case J ■ 2m + 1, 
are all real positive. 

The first result could be derived directly for the block matrix with 
the same structure. The second one Implies that a necessary condition for 
the matrix Ij to have only real roots. Is that the coefficients of P„(x), 




I- 

i- 


r ' 








no 


in case J “ 2m, or Q^(x)» case J - 2m + 1, alternate In sign, or 
equivalently, that the coefflclente of 0 *“’ or 6„ \ reapectlvely, glve.1 by 


„(") . (-l)"-''pi'')(0)/vl ' 

► 

where v ■ 0, 1, 2, . . . , m, be all positive. Recurrence formulas 
these coefficients can easily be obtained using Equation (Gil). In 
lar, setting x ™ 0 in this equation and multiplying the result by 
gives : 


(G17) 


for 

partlcu- 


(- 1 ) 


mfl 


(o) (o) 

uri-l anri"! 2mf2 m 


^mfl “nrfl ®2nH-2®2iiri-3^iii 


(G18a) 

(G18b) 


Similarly, dlfferenclatlng Equation (Gll) v times (v > 1) with respect to 


nrt-l-v , 


X, setting x ■ 0, and multiplying the result by (-1) /v! give: 

(v) (v-1) (v) 

Vl-^m + "2mfl^2»f2“« 


„(v) . (v) 

‘^nrfl mfl 2m4-2 2m+3'^m 


(V) 


(G19a) 

(G19b) 


which in fact, in view of Equation (G18), can be applied with v - 0 if one 

defines 8^”^^ = 0. Equation (G19) should be completed by the following 
m 

"boundary" conditions: 

_ 1 (G20) 

m IQ 

which simply state that x® is the leading term of both 

From this, it is easy to derive explicit formulas for ®nd 

for V - 0 and .1, now denoted more simply by a , 8 , a' , and 6', respec- 
tlvely. In particular, applying Equation (G18a) recursively gives: 


< 


1 




! i 


^ ^ i ) 


.■. ! 


, , - , , -. . . (.’■•■ ■» --1 ' -r- i -*-i' 
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m 2m- 1 2m m-l 


A A A A Ot 

2m- 3 2m-2 2m- 1 2m m-2 




- . . 


. a„ a 
2m o 


2m 


(G21) 


where Eqtiatlon (G20) has been used with m > 0. Now, write Equation (G18b) 
with m - k - 1 (k ■ 1, 2, . . . , m) and multiply the result by 


®2k+2®2k+3 


. a„ . , . This gives 
2mi-l ® 


®rk+2®2k+3 * ' * ®2mfl\ " *2k“2k+l * ' ' ®2mfl\-l 


A . A. * « • A. , • 

1 2 2m4-l 


a 


2k+l 


(G22) 

where Equation (G21) has been used. (For k ■ m, the coefficient of 
in this equation is understood to he one.) Then, write Equation (G22) for 
k - 1, 2, . . . , m, add the resu.’.ting relationships, use Equation (G20) with 
m » 0, and simplify the result to get: 

m 

> 1 

(G23) 


“ a, a. . 

m 12 


2mrfl 


^ a.,. 


k"0 


2k+l 


Similarly, write Equation (G19a) with m ■ k - 1 (k ■ 1, 2, . . . , m) and 


V ■ 1, and multiply the result by ®2k+l*2k+2 


. a^^ to get: 


A 


2k+l 2k+2 


• • • A ^ 


a_a’ . + 


®1®2 • 


2m 


2m k 2k-l 2k 2m k-1 


2k 


ic 

E 


‘2j-l 

(G24) 


.iT. »-"f ' rr- 
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where Equation 
k - 1, 2, . . . 


(G23) has been used. Then, write Equation (G24) for 
, m and add the resulting relationships to get: 


m k 



(G25) 


where the fact that a’ 
with (1*1, 

®2i+2®2£+3 • ' • ®2 qH 


-■ 0 has been used. Finally, write Equation (G19b) 
2, . . ,, m) and v >■ 1, and multiply the result by 
to get: 


^2 £+2^2 £+3 * 


®2nrfl®£ “ ®2£^2£+l ' ’ ' ®2mfl®£-l 


+ 



where Equation (G25) has been used. Then write Equation (G26) for 
£ ■• 1, 2, . . m and add the resulting relationships to get: 

m £ k 


6 


t 

m 





—T—T — 


(G26) 


(G27) 


where the fact Chat » 0 has been used. 

In view of the Equations (G21), (G23), (G25), and (G27), it appears 
that in the event one or several alternations of sign occur In the sequence 
Sj, without separation of positive from negative values by at least one zero. 
It Is most unlikely that the coefficients or are all positive. 

It suffices that only one of these coefficients be negative, for ?„(*) or 
Q„,(x) to have at least one root x* which is not real positive; then, the 
matrix Fj has, in particular, the eigenvalues these are complex and 

one of the two has a positive Imaginary part; to that one corresponds an 
eigenvalue of the loatrix (sae Equation (70)) with a negative real part 
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which has a destabilizing effect, sii it acts like a negative smoothing. 
More precisely, suppose for example, that J ■ 2m and that “n 0* Thi3 
occurs for example when a^, s. 2 , • • •, ^2r negative for some odd 

value of r. Then let r^, . . • be the real roots of Pj,(x) and 

zj, ij, zj, Z 2 , ... be its complex roots. The coefficient is 

siti^>ly the product of these roots so that: 


“m 


rirj . 


. (Zj5j)(z2Z2) . . 


- rir2 . 


. . . < 0 


(G28) 


This shows that Pjg(x) then has an odd nuafl>er of real negative roots 
n^, n 2 > . . . to which correspond an odd number of pairs of purely imaginary 
eigenvalues l/-nj and -iZ-n^ , l/^ 2 » -l»^-n 2 , • . . for the matrix F 

and the real eigenvalues -/-nj and /-nj*, -/-nj and /-n 2 , . . . for the 


matrix C., 


Remark 1: 


A similar situation occurs if J ■ 2m + 1 and 8^, < 0, 


The fact that the eigenvalues of tlie matrix Fj are real for Cases I 
and II can be derived from Equation (GIO) and (Gll). In fact, this result was 
originally obtained in this way. In that first analysis, the following 
separation properties were found to be true for Case I: 

0 < xj « yj < xj < y 2 < . . . < x„ < y„ 

0 < xj < Xj < x^ < X 2 < . . . < x; < x„ < x^j J (G29) 

0 < y| < y < y’ < y < . . . < y' < v < v* 

■^1 ■’^2 ^2 ^m 

where {x^!, {y^} (j - 1, 2, . . m), {xj} and {yj} (J - 1, 2, . . . , m + 1) 

are the roots of P„(x), Q„(*:), P^j (x) and <^j(x), respectively. The 
derivation of Equation (G29) is omitted here, since no r>artlcular applica- 
tion of it was found. 
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However, Equation (G29) shows that the eigenvalues of the matrix Tj 
and thus of the matrix are simple, which Implies that both matrices 

are dlagonallzable, an already known fact. 

Remark 2 ; 

Consider the case where the Euler implicit method Is applied to the 
following one-dlmenslonal (generalized) wave-equation: 

+ Ia(x,t)u]^ - 0 (G30) 

If no smoothing is applied, the solution-vector u”"*”^ at the n + 1st 
time-step is given by: 


where the following definitions are made: 

Cf ' . irld 0, .Jti) 

- iDr”"'^D 
D - Dlag(i^) 


(G31) 


(G32) 


Suppose that a uniform bound on u^ is required. For this, assume that 
a(x,t) >0 so that the results of Case I are applicable. In particular, 
the transformations that diagonalize the matrices f and c""*"^ (or 
L) are, respectively, S, VG, and X ■ D7G where G Is some orthogonal 
matrix, and the diagonal siatrlces 7 and D are given by Equations (C4) 
and (G32), respectively. Consider first the case where a(x,t) does not 
depend on time, so that L does not either. Then 

lu"! - l(XAX- jVl - IIxaV1u°I 


• t • -rf ■■ "1 • .'-"t-— -(•'-'-''•■fi'- ' 
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» 


f 


where v * Ixl • Hx"^l is a condition number, and A Is the diagonalized 
form of the matrix L whose eigenvalues are given by; 


X 


J 


I + lYj 


(G34) 


where Yj 
that I >,^ I 


is an eigenvalue of the matrix Fj, that is, a real number, so 
< 1. If the euclidean norm is selected, lAl < 1 so that 

lu”ll < vlu°l (G35) 


which is the bound we were looking for. The value of v has no Importance 
In this case. 

However, if now a(x,t) does depend on time, it appears inevitable to 

use the following bound for L: 

iLl - IXAX“1| < V (G36) 

where v is an upper bound for v (which depends on n). This gives; 

lu”l < v"lu°l (G37) 

It hence appears of interest to evaluate v. For this, recall that fl Is 

orthogonal and D unitary, so that: 

IIXl « iDVslII = II Vl - Max|6 

J ^ 

lx~l| « lfi*^7"li)ll - llV'^n - 1/Mln|6 "I 

j ^ 

where euclidean norm has been used. This gives 



V * Sup[Max| 6 /Mln|fi ”1 ] 
n j J J ^ 

- SupJMjixla "|/Mln|a "I (G39) 

n f :i ^ J ^ 

where Equation (G3) has been used. Clearly, v > 1 unless a(x,t) only 
depends on time, and no uniform bound for u*' can be derived from Equa- 
tion (G37), which, to be rigorous, only reveals the failure of this attempt. 





I 




V ! i 

4 . . *- 
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,4 


r 
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However, for a practical problem, one anticipates v to be very large If 
the eigenvalues a^^ are themselvos subject to large variations, :his 
suggests that despite the fact that Equation (G36) Is a conservative esti- 
mate, whe operator L of Equation (G31) is unlikely to be contracting as 
one would wish in order to apply the contraction-mapping theorem cited in 
Section IIA (i,u., contracting in the same norm sense for all n), 
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